I think I got a few simple concepts confused although the code below will
eventually do what I need, I’m not sure it is the best way to do.
Any comments? Thanks.
load "isoline"
load "medit"
// What does "Curve()" return and what is "P" in boarder?
// What type of object does boarder actually return ?
// I wanted to add parameters for location to a box building
// boarder but could not find a good way to do that.
// You can see the AntBox solution below.
// I could not come up with a macro I could pass to buildmesh
// and could just assign the boarder to some new object.
// How can I define a surface for integration ( I ultimaately
// want the empty boxes to be antenna elements and will integate
// the normal Poynting vectorS \cdot \hat{n} over a bounding
// surface ) without having
// it mess up the mesh? I guess "intalledges" does not do that :)
// as currently written, the mesh this code display is about
// the required result but the radiation pattern circle
// is includes in the meshing. Also I wanted to plot
// normal poynting along the contour.
real[int,int] xy(3, 1); //to store the isoline points
int[int] be(1); //to store the begin, end couple of lines
real[int,int] xye(3, 1); //to store the isoline points
int[int] bee(1); //to store the begin, end couple of lines
real Lx=1;
real Ly=1;
func int pcurce(real[int,int] xy, int[int] be, int nbc, real[int] viso )
{
for (int c = 0; c < nbc; ++c){
int i0 = be[2*c], i1 = be[2*c+1]-1;
cout << "Curve " << c << endl;
for(int i = i0; i <= i1; ++i)
cout << "x= " << xy(0,i) << " y= " << xy(1,i) << " s= " << xy(2, i) << endl;
plot([xy(0, i0:i1), xy(1, i0:i1)], wait=true, viso=viso, cmm=" curve "+c);
}
return 0;
} // pcurce
{
func real prect(real xx, real yy, real asp) {
real n=50;
real un=1.0/n;
//return pow((xx*asp/Lx)^n+(yy/Ly)^n, un);
real xa=(xx<0)?(-xx):xx;
real ya=(yy<0)?(-yy):yy;
if ((xa*asp)>(ya) ) return (xa*asp); else return (ya);
}
mesh Th = square(10, 10);
mesh Thc = square(10, 10,[Lx*(x-.5),Ly*(y-.5)]);
fespace Vh(Th, P1);
fespace Vhc(Thc, P1);
Vh u = sqrt(square(x-0.5) + square(y-0.5));
Vhc w= prect(x,y,.5);
real iso = 0.2 ;
real[int] viso = [iso];
//plot(u, viso=viso,Th);//to see the iso line
//plot(w, viso=viso,fill=1,wait=1);//to see the iso line
//medit("asdf",Thc,w);
cout<<w<<endl;
int nbc = isoline(Th, u, xy, close=1, iso=iso, beginend=be, smoothing=0.1);
int nbce = isoline(Thc, w, xye, close=1, iso=iso, beginend=bee, smoothing=0.1);
cout << "Number of the line component = " << nbc << endl;
cout << "Number of points = " << xy.m << endl;
cout << "be = " << be << endl;
//int xxx= pcurce( xye, bee, nbce,viso);
// shows the lines component
if (false) {
for (int c = 0; c < nbc; ++c){
int i0 = be[2*c], i1 = be[2*c+1]-1;
cout << "Curve " << c << endl;
for(int i = i0; i <= i1; ++i)
cout << "x= " << xy(0,i) << " y= " << xy(1,i) << " s= " << xy(2, i) << endl;
plot([xy(0, i0:i1), xy(1, i0:i1)], wait=true, viso=viso, cmm=" curve "+c);
}
} // false
}
cout << "length of last curve = " << xy(2, xy.m-1) << endl;
border Curve0(t=0, 1){
int c=0; //component 0
int i0=be[2*c], i1=be[2*c+1]-1;
P=Curve(xy, i0, i1, t); //Curve 0
label=1;
}
border Curve1(t=0, 1){
int c=1; //component 1
int i0=be[2*c], i1=be[2*c+1]-1;
P=Curve(xy, i0, i1, t); //Curve 1
label=1;
}
border Curve3(t=0, 1){
int c=0; //component 0
int i0=bee[2*c], i1=bee[2*c+1]-1;
P=Curve(xye, i0, i1, t); //Curve 0
label=3;
}
border Curve4(t=0, 1){
int c=1; //component 1
int i0=bee[2*c], i1=bee[2*c+1]-1;
P=Curve(xye, i0, i1, t); //Curve 1
label=3;
}
//real ytop=.1; real ybot=0; real xleft=-.8; real xright=.8;
// cout<<"t="<<t<<endl;
macro AntBox(CurveName,l,r,top,b,lab)
border CurveName(t=0, 1){
real ytop=top; real ybot=b; real xleft=l; real xright=r;
if (t<.45)
{ y=ybot; x=(xright-xleft)/.45*t+xleft;; }
else if (t<.5) { x=xright; y=-(ybot-ytop)/.05*(t-.45)+ybot; }
else if (t<.95) { y=ytop; x=xright-(xright-xleft)/.45*(t-.5); }
else { x=xleft; y=(ybot-ytop)/.05*(t-.95)+ytop; }
;label=lab;
} // EOM
// Curve5
//cout<<"x "<<x<<" y= "<<y<<" t="<<t<<endl;
//macro LPDA(t,b,l,r,n)
//mesh Th = buildmesh(Curve3(100)+Curve4(100));
//mesh Th = buildmesh(Curve5(100)+Curve4(100));
//border x1=Curve5(100);
//mesh Th = buildmesh(Curve5(100));
AntBox(Curve5,-.8,.8,.1,0,100);
AntBox(Curve7,-.6,.6,-.3,-.4,101);
AntBox(Curve6,-2,2,1,-1,102);
int cl=99;
border rads(t=0, 2.*pi){x=.9*cos(t); y=.9*sin(t); label=cl; }; // label=cl };
// seg faults
//mesh Th =emptymesh(buildmesh(rads(100)))+ buildmesh(Curve6(100)+Curve5(-100)+Curve7(-100));
//mesh Th =buildmesh(Curve6(100)+Curve5(-100)+Curve7(-100));
mesh Th =buildmesh(Curve6(100)+rads(20)+Curve7(-100) );
//mesh Th = Curve5(-100)+buildmesh(Curve6(100)+Curve7(-100));
//Th=change(Th,fregion=((x*x+y*y)<=.25)?99:region);
fespace Vh(Th,P2);
Vh hh=x*x;
//real v=int1d(Th,cl)(hh*N.x);
//real v=int1d(Th,cl)(1.0);
//real v=intalledges(Th,99)(hh*N.x);
//real v=intalledges(Th,99)(1.0);
real v=int1d(Th,99)(1.0)/.9;
cout<<" v = "<<v<<endl;;
//ytop=1; ybot=-1; xleft=-1; xright=1;
//Th = Th + buildmesh(Curve5(100));
plot(Th, wait=true);
if (false) {
plot(Curve1(100)); //show curve
mesh Th = buildmesh(Curve1(-100));
plot(Th, wait=true);
Th = buildmesh(Curve0(100)+Curve1(100));
plot(Th, wait=true);
} // false
[isoline.edp|attachment](upload://tHrajxFjRydMwzT6SRh2nvceAq3.edp) (5.0 KB)