Die elementaren Grundlagen zum Integralsatz werden in den Videos unten genannt, ich setze sie als bekannt vorraus...
Nun zu deinem speziellen Problem:
Du solltest dich irgendwann im Studium mal mit den sogenannten "Normen" beschäftigt haben, und hoffentlich hat man dir auch dieses bzw. ein ähnliches Bild gezeigt:

Hier sind in GeoGebra folgende Normen (bzw. nenne ich sie in Analogie so - unabhängig davon, ob sie die Definition exakt erfüllen) für r=4 dargestellt (von außen nach innen):
Unendlichkeitsnorm (das Quadrat) mit Fläche `A=4*r^2`, 6-Norm, 3-Norm, 2-Norm, also der Kreis mit der bekannten Gleichung: `x^2+y^2=r^2` bzw. `y=+-(r^2-x^2)^(1/2)` und Fläche `A=pi*r^2`, 4/3-Norm, 1-Norm also `abs(x)+abs(y)=r` bzw. `y=+-(r-abs(x))` (gekipptes kleines Quadrat) mit der Fläche `A=2*r^2`, 2/3-Norm (dein Fall - rot markiert), 1/2-Norm, 1/3-Norm, 1/4-Norm - wobei jeweils nur der Teil innerhalb des großen Quadrates zu betrachten ist!

Diese Dinge begegnen einem immer mal wieder sporadisch und man sollte sie im Hinterkopf behalten...
Das kleiner gleich Zeichen beschreibt jetzt den gesamten inneren Bereich der Kurve:
`y=+-sqrt(1-3x^((2)/(3))+3x^((4)/(3))-x^(2))` Wäre die Gleichung deiner Randkurve für r=1 oder allgemein:
`y=+-sqrt(r^2-3*r^(4/3)*x^((2)/(3))+3*r^(2/3)*x^((4)/(3))-x^(2))`
Ich habe hier die binomische Formel für hoch 3 angewandt, um dir den `r^2-x^2` Teil sichtbar zu machen!

Jetzt können wir natürlich auch noch eine allgemeine Funktion angeben:
`y=+-(r^(p/q)-abs(x)^(p/q))^(q/p)`
Diese einfache Gleichung beschreibt die Randkurve jeder positiv rationalen Norm, sofern man den Bereich von `x` auf `[-r,r]` einschränkt. Mit den richtigen Fallunterscheidungen lassen sich auch die Beträge auflösen. Es ist nun so, dass leicht eingesehen werden kann, dass jede dieser Funktionen symmertrisch zu beiden Achsen des Koordinatensystem ist, und sich die durch die Randkurve begrenzte Fläche aus 4 mal der Fläche im ersten Quadranten zusammensetzt.
Das natürliche Vektorfeld des Koordinatensystems ist nun einfach `V=((x/2),(y/2))` - Mit der Divergenz `1/2+1/2=1`. Für eine Flächenberechnung integrieren wir ja stets über eine (gedachte) 1, also z.B. `A=intdA*1=intint"div"(vec(v))dxdy`.

Auch wenn du es nicht musst - zeigen wir nun, dass folgende Darstellungen äquivalent sind (zunächst nur für den ersten Quadranten also mit `tin[0,pi/2]` und `sin(x)=abs(sin(x))` sowie `cos(x)=abs(cos(x))` , daher ohne Beträge):

I) `abs(x)^(p/q)+abs(y)^(p/q)=r^(p/q)`
II) `y=+-(r^(p/q)-abs(x)^(p/q))^(q/p)` --> Folgt aus trivialem Umstellen von I)
III) `(r*cos(t)^((q*2)/p),r*sin(t)^((q*2)/p),t,0,2pi)`

Setzen wir III) in I) ein, so erhalten wir:
`(r*cos(t)^((q*2)/p))^(p/q)+(r*sin(t)^((q*2)/p))^(p/q)=r^(p/q)`
`r^(p/q)*((cos(t)^((q*2)/p))^(p/q)+(r*sin(t)^((q*2)/p))^(p/q))=r^(p/q)`
`r^(p/q)*(cos(t)^2+sin(t)^2)=r^(p/q)`

Dies ist offensichtlich eine wahre Aussage.
Weitere Untersuchungen führen darauf, dass für `p=2` sowie `q` ungerade und `t` im Bereich von `[0,2pi]` tatsächlich nicht nur der erste Quadrant, sondern die gesamte Figur dargestellt wird. Somit können wir auch die Parametrisierung guten Gewissens verwenden.

Darstellung des Randes deiner Fläche für `r=4`

Überlegen wir uns nun, wie der Normalenvektor `vecv_n` aussieht:
Dazu bestimmen wir die Ableitungen nach t:
x: `d/(dt)r*cos(t)^3=r*3*cos(t)^2*(-sin(t))`
y: `d/(dt)r*sin(t)^3=r*3*sin(t)^2*cos(t)`
Nun vertauschen wir x und y und ändern das Vorzeichen von x:
`vec(v_n)=((r*3*sin(t)^2*cos(t)),(r*3*cos(t)^2*sin(t)))`

Also Nächstes bestimmen wir den Betrag (mit mehr Erfahrung könnten wir diesen Schritt auch weglassen, da sich der Betrag am Ende wieder wegkürzen wird):
`abs(vec(v_n))=sqrt(r^2*9*(sin(t)^4*cos(t)^2+cos(t)^4*sin(t)^2)`
`=3r*sqrt((sin(t)^2*cos(t)^2*(sin(t)^2+cos(t)^2)))`
`=3r*sqrt((sin(t)^2*cos(t)^2))`
`=3r*abs(sin(t)*cos(t))`

Also erhalten wir:
`vec(v_n)/abs(vec(v_n))=1/(3r*abs(sin(t)*cos(t)))*((r*3*sin(t)^2*cos(t)),(r*3*cos(t)^2*sin(t)))`
`=((sin(t)),(cos(t)))` --> Erster Quadrant und Dritter Quadrant
`=((-sin(t)),(-cos(t)))` --> Zweiter Quadrant und Vierter Quadrant
Mit etwas Erfahrung hätten wir dies auch gleich einsehen können.

Hier nochmal als Bildlich Darstellung - Normalenvektor auf der Randkurve und dargestellt durch Vektoren auf dem Einheitskreis.

Wir müssen nun folgendes Integral berechnen (unter Ausnutzung der Symmetrie "4 mal der erste Qudadrant", also nur bis `(2pi)/4=pi/2` integrieren, um das Ganze etwas abzukürzen):

`A=4*int_0^(pi/2)ds((x/2),(y/2))*vec(v_n)/abs(vec(v_n))` | Jetzt wie immer "`ds=abs((ds)/(dt))*dt`" verwenden...
Hier kürzt sich auch die Länge wieder heraus, weil `abs((ds)/(dt))=abs(vec(v_n))`
`=4*int_0^(pi/2)dt((x/2),(y/2))*((r*3*sin(t)^2*cos(t)),(r*3*cos(t)^2*sin(t)))` |x und y durch unsere Parametrisierung ersetzen!
`=4*int_0^(pi/2)dt(r/2*cos^3(t)*(r*3*sin(t)^2*cos(t))+r/2*sin^3(t)*(r*3*cos(t)^2*sin(t)))`
`=4*int_0^(pi/2)dt(r^2*3/2*(cos^4(t)*sin^2(t)+sin^4(t)*cos^2(t)^2))`
`=4*int_0^(pi/2)dt(r^2*3/2*(cos^2(t)*sin^2(t)*(sin^2(t)+cos^2(t))))`
`=4*int_0^(pi/2)dt(r^2*3/2*cos^2(t)*sin^2(t))`

Dieses Integral ist nun zu lösen:
Es gelten die bekannten Identitäten:
`cos^2(x)=1/2*(1+cos(2x))` und `sin^2(x)=1/2*(1-cos(2x))`
Und somit:
`sin^2(x)*cos^2(x)=1/4*(1-cos(2x)^2)=1/4*(1-1/2*(1+cos(4x)))`
`=1/4-1/8-1/8*cos(4x)=1/8-1/8*cos(4x)`

`A=4*int_0^(pi/2)dt*(3/2*r^2)*(1/8-1/8*cos(4x))`
`=4*3/2*r^2*[1/8*x-1/32*sin(4x)]_0^(pi/2)` | mit `sin(4*pi/2)=sin(2pi)=sin(0)=0`
`=4*3/2*r^2*[(pi/2)/8]=(4*3*pi*r^2)/32=3/8*pi*r^2`
 
Für `r=4` erhalten wir z.B.: `A=6*pi~~18.849556`
Eine numerische Integration liefert für diesen Fall: `A~~18.849556`
Wir können also von der Richtigkeit des Ergebnisses ausgehen!