Hello, the code is always stopping at the prob command and i think the problem comes from the mesh creation of the spheres but still can not detect where.
Here is the code i am working on:
load “msh3” // provides functionality for reading and manipulating mesh files in the .mesh format, which is commonly used for storing 2D and 3D meshes. By loading this module, you gain access to functions and commands related to mesh manipulation and processing.
load “medit” // loads the visualization module named “medit”. The “medit” module in FreeFEM is used for visualizing finite element meshes and their associated data. It allows you to visualize the mesh geometry, boundary conditions, solution fields, and other relevant information in a graphical interface.
load “tetgen”
searchMethod=1; //sets the value of the variable searchMethod to 1
verbosity=1; //sets the value of the variable verbosity to 1
// Defining the elastic constants
real c111=6530, c121=3870, c131=3870, c141=0, c151=0,c161=0;//MPa
real c211=3870, c221=6530, c231=3870, c241=0,c251=0,c261=0;
real c311=3870, c321=3870, c331=6530, c341=0,c351=0,c361=0;
real c411=0, c421=0, c431=0, c441=1330,c451=0,c461=0;
real c511=0, c521=0, c531=0, c541=0,c551=1330,c561=0;
real c611=0, c621=0, c631=0, c641=0,c651=0,c661=1330;
//material of fiber
// Defining Elastic Constants
real c112=6530, c122=3870, c132=3870, c142=0, c152=0,c162=0;//MPa
real c212=3870, c222=6530, c232=3870, c242=0,c252=0,c262=0;
real c312=3870, c322=3870, c332=6530, c342=0,c352=0,c362=0;
real c412=0, c422=0, c432=0, c442=1330,c452=0,c462=0;
real c512=0, c522=0, c532=0, c542=0,c552=1330,c562=0;
real c612=0, c622=0, c632=0, c642=0,c652=0,c662=1330;
real c113=755430, c123=181490, c133=255770, c143=0, c153=0,c163=0;//MPa
real c213=181490, c223=755430, c233=255770, c243=0,c253=0,c263=0;
real c313=255770, c323=255770, c333=1641900, c343=0,c353=0,c363=0;
real c413=0, c423=0, c433=0, c443=1222680,c453=0,c463=0;
real c513=0, c523=0, c533=0, c543=0,c553=1222680,c563=0;
real c613=0, c623=0, c633=0, c643=0,c653=0,c663=286190;
// Declaring and initizalizing variables (young Modulus, poisson,…)
real test=116800,b=6290;
real E2in,nuin;
real wtot,w=3139.32,rint;
real R = 1;
real R1 = R + 0.5;
int nz=3;
while(1){ // an infinite loop as long as the condition is true
c112=6530; c122=3870;c132=3870; c142=0; c152=0;c162=0;//MPa
c212=3870; c222=6530; c232=3870; c242=0;c252=0;c262=0;
c312=3870; c322=3870; c332=6530;c342=0;c352=0;c362=0;
c412=0; c422=0; c432=0; c442=1330;c452=0;c462=0;
c512=0; c522=0; c532=0; c542=0;c552=1330;c562=0;
c612=0; c622=0; c632=0; c642=0;c652=0;c662=1330;
int left=1,right=2,front=3,back=4,l1=5,l0=6;
real a=3, d=0.5, h=3;
border b1(t=0.5,-0.5) {x=at; y=-a/2; label=back;}; //defines a border named b1. The border is parametrized by t with limits from 0.5 to -0.5. The x and y coordinates of points on this border are defined as functions of t
border b2(t=0.5,-0.5) {x=a/2; y=at; label=right;}; // label allows us to give the border a description or an identifier to refer to it later
border b3(t=0.5,-0.5) {x=at; y=a/2; label=front;};
border b4(t=0.5,-0.5) {x=-a/2; y=at; label=left;};
real cx = 0;
real cy = 0;
real cz = 1.5;
mesh Th=square(10,20,[xpi-pi/2,2y*pi]);
// a parametrization of a sphere
func f1 = cx + cos(x)*cos(y);
func f2 = cy+ cos(x)*sin(y);
func f3 = cz+ sin(x);
real czo = 1;
func f4 = cx + cos(x)*cos(y);
func f5 = cy+ cos(x)*sin(y);
func f6 = czo+ sin(x);
func f1x=sin(x)cos(y);
func f1y=-cos(x)sin(y);
func f2x=-sin(x)sin(y);
func f2y=cos(x)cos(y);
func f3x=cos(x);
func f3y=0;
// M = DF^t DF
func m11=f1x^2+f2x^2+f3x^2;
func m21=f1xf1y+f2xf2y+f3xf3y;
func m22=f1y^2+f2y^2+f3y^2;
verbosity=10;
real[int] domain =[0.,0.,0.,1,0.01];
mesh3 Th3=tetgtransfo(Th,transfo=[Rf1,Rf2,Rf3],nbofregions=1,regionlist=domain);
mesh3 Th4=tetgtransfo(Th,transfo=[R1f4,R1f5,R1*f6],nbofregions=1,regionlist=domain);
int reg0=3 ,reg1=2,reg2=1;
int nnb=7, nni=10;
mesh Thx=buildmesh(b1(-nnb)+b3(nnb)+b2(-nnb)+b4(nnb));// specifies the boundaries and interfaces used to build the mesh.
// nnb,20 are variables that controls the number of points or subdivisions along the boundary curves
int[int] REG=[Thx(0,a/0.49).region,reg0, Thx(0,0.5).region,reg1,Thx(0,0).region,reg2];
// REG is an integer array used to define regions within the mesh Thx
// Thx(0,a/0.49).region: This part of the code extracts the region number of the mesh Thx at a specific point. The point is defined by the coordinates (0, a/0.49) in the domain.
cout << reg0 << " " << reg1 << " " << reg2 << endl; // prints the values of the integer variables reg0, reg1, and reg2 to the standard output stream
// end1 ensures that each set of values is printed in a new line
//int[int] REG=[reg0,10,reg1,11,reg2,12];
Thx=change(Thx,region=REG); // The change function is used to modify the mesh Thx by assigning regions to its elements based on the values specified in the REG array.
//cout << regions(Thx)<< endl;
plot(Thx,wait=1); ;
{ // for cleanning memory…
int[int] old2new(0:Thx.nv-1); //This creates an integer array old2new with indices ranging from 0 to the number of vertices (Thx.nv) minus one. This array will map old vertex indices to new vertex indices.
fespace Vh2(Thx,P1); //This defines a finite element space Vh2 on the mesh Thx using P1 elements
Vh2 sorder=x+y; //Here, sorder is assigned the sum of the x and y coordinates of the vertices. This will be used to determine the sorting order of the vertices.
sort(sorder,old2new); //The sort function sorts the vertices based on sorder and fills old2new with the new indices such that the vertices are sorted in ascending order of x + y
int[int] new2old=old2new^-1; // inverse the permuation
//for(int i=0;i< Th.nv;++i) // so by hand.
// new2old[old2new[i]]=i;
Thx= change(Thx,renumv=new2old); //The change function is used to apply the new vertex numbering to the mesh Thx, using new2old to map old vertex indices to new ones.
sorder=0:Thx.nv-1; //This reinitializes sorder with values from 0 to the number of vertices minus one, effectively resetting it.
}
{
fespace Vh2(Thx,P1);
Vh2 nu;
nu=0:Thx.nv-1;
plot(nu,cmm=“nu=”,wait=1);
}
mesh ThHs = Th + Thx;
real haut=10;
real bas=11;
int[int] rup=[reg0,haut,reg1,haut,reg2,haut], rlow=[reg0,bas,reg1,bas,reg2,bas], rmid=[0,6], rtet=[0,41];
// Arrays rup, rlow, rmid, and rtet are defined. These arrays seem to store pairs of region identifiers and their associated values (e.g., haut and bas).
func zmin=0;
func zmax=h;
mesh3 mohamad=buildlayers(Thx, nz, zbound=[zmin,zmax], // A 3D mesh named mohamad is created from a 2D mesh Thx with nz layers between the z-boundaries defined by zmin and zmax. The mesh faces are referenced using rup for the top and rlow for the bottom.
reffaceup=rup, reffacelow=rlow);
for(int i=1;i<=6;++i) // loop calculates and prints integrals over the mesh faces for regions 1 to 6.
cout << " int " << i << " : " << int2d(mohamad,i)(1.) << " " << int2d(mohamad,i)(1./area) << endl;
savemesh(mohamad,“Th3.mesh”); // The mesh mohamad is saved to a file named Th3.mesh
mesh3 Thmhmd = Th3 + Th4;
plot(Thmhmd,wait=1);
mesh3 Th5 = Thmhmd + mohamad;
plot (Th5, wait=1);
fespace Wh(Th5,P1); //Wh is defined as a finite element space on the 3D mesh mohamad using linear elements (P1).
varf a1(u,v)=on(haut,u=1); //Two variational forms, a1 and a2, are defined for imposing Dirichlet boundary conditions on the faces identified by haut and bas
varf a2(u,v)=on(bas,u=1);
Wh u=0; //The function u is defined in the finite element space Wh and initialized to 0.
u=a1(0,Wh,tgv=1); //variational problem is solved with a1, setting u to 1 on the haut boundary. The result is stored in the array u, and the function u is plotted with the label “label 1”.
plot(u,wait=1,cmm=“label 1”);
u=a2(0,Wh,tgv=1);
plot(u,wait=1,cmm=“label 3”); //second variational problem is solved with a2, setting u to 1 on the bas boundary. The new result is stored in u, and u is plotted again with the label “label 3”.
fespace Vh66(Th5,[P2,P2,P2], periodic=[ [back,x,z],[front,x,z] ,[left,y,z],[right,y,z] ] );
// defines a finite element space Vh66 on the 3D mesh mohamad using quadratic elements (P2) for each spatial dimension (x, y, z). The periodic keyword specifies the periodic boundary conditions.
fespace fspace0(Th5,P0);
while(1)
{
// We define the coefficient arrays c111, c112, etc., and a function initializeCoefficients() to initialize them.
// The function getCoefficient(i, j, regionId) returns the appropriate coefficient array based on i, j, and regionId
fspace0 c11=(region==reg0)*c111+(region==reg1)c112+c113(region==reg2);
fspace0 c12=(region==reg0)*c121+(region==reg1)c122+c123(region==reg2);
fspace0 c13=(region==reg0)*c131+(region==reg1)c132+c133(region==reg2);
fspace0 c14=(region==reg0)*c141+(region==reg1)c142+c142(region==reg2);
fspace0 c15=(region==reg0)*c151+(region==reg1)c152+c153(region==reg2);
fspace0 c16=(region==reg0)*c161+(region==reg1)c162+c163(region==reg2);
fspace0 c21=(region==reg0)*c211+(region==reg1)c212+c213(region==reg2);
fspace0 c22=(region==reg0)*c221+(region==reg1)c222+c223(region==reg2);
fspace0 c23=(region==reg0)*c231+(region==reg1)c232+c233(region==reg2);
fspace0 c24=(region==reg0)*c241+(region==reg1)c242+c243(region==reg2);
fspace0 c25=(region==reg0)*c251+(region==reg1)c252+c253(region==reg2);
fspace0 c26=(region==reg0)*c261+(region==reg1)c262+c263(region==reg2);
fspace0 c31=(region==reg0)*c311+(region==reg1)c312+c313(region==reg2);
fspace0 c32=(region==reg0)*c321+(region==reg1)c322+c323(region==reg2);
fspace0 c33=(region==reg0)*c331+(region==reg1)c332+c333(region==reg2);
fspace0 c34=(region==reg0)*c341+(region==reg1)c342+c343(region==reg2);
fspace0 c35=(region==reg0)*c351+(region==reg1)c352+c353(region==reg2);
fspace0 c36=(region==reg0)*c361+(region==reg1)c362+c363(region==reg2);
fspace0 c41=(region==reg0)*c411+(region==reg1)c412+c413(region==reg2);
fspace0 c42=(region==reg0)*c421+(region==reg1)c422+c423(region==reg2);
fspace0 c43=(region==reg0)*c431+(region==reg1)c432+c433(region==reg2);
fspace0 c44=(region==reg0)*c441+(region==reg1)c442+c443(region==reg2);
fspace0 c45=(region==reg0)*c451+(region==reg1)c452+c453(region==reg2);
fspace0 c46=(region==reg0)*c461+(region==reg1)c462+c463(region==reg2);
fspace0 c51=(region==reg0)*c511+(region==reg1)c512+c513(region==reg2);
fspace0 c52=(region==reg0)*c521+(region==reg1)c522+c523(region==reg2);
fspace0 c53=(region==reg0)*c531+(region==reg1)c532+c533(region==reg2);
fspace0 c54=(region==reg0)*c541+(region==reg1)c542+c543(region==reg2);
fspace0 c55=(region==reg0)*c551+(region==reg1)c552+c553(region==reg2);
fspace0 c56=(region==reg0)*c561+(region==reg1)c562+c563(region==reg2);
fspace0 c61=(region==reg0)*c611+(region==reg1)c612+c613(region==reg2);
fspace0 c62=(region==reg0)*c621+(region==reg1)c622+c623(region==reg2);
fspace0 c63=(region==reg0)*c631+(region==reg1)c632+c633(region==reg2);
fspace0 c64=(region==reg0)*c641+(region==reg1)c642+c643(region==reg2);
fspace0 c65=(region==reg0)*c651+(region==reg1)c652+c653(region==reg2);
fspace0 c66=(region==reg0)*c661+(region==reg1)c662+c663(region==reg2);
real E11,E22,E33,E23,E13,E12,p1,p2,p3; //p1 and p2 are macro electric fields
real k111,k112,k113,k221,k222,k223,k331,k332,k333,k231,k232,k233,k131,k132,k133,k121,k122,k123;
real S111,S112,S113,S221,S222,S223,S331,S332,S333,S231,S232,S233,S131,S132,S133,S121,S122,S123;
real G11,G22,G33,G23,G13,G12;
macro u [ux,uy,uz] //
// macro defines u to be a shorthand for the list [ux, uy, uz].
macro v [vx,vy,vz] //
Vh66 u,v; // specifies that u and v are vector fields in the finite element space Vh66. This declaration can be used to define and solve variational problems involving vector fields.
fspace0 epsilonx,epsilony,epsilonxy,epsilonz,epsilonxz,epsilonyz,sigmax,sigmay,sigmaxy,sigmaz,sigmaxz,sigmayz,px,py,pz;//
func C6x6 = [
[c11,c12,c13,c14,c15,c16 ],
[c21,c22,c23,c24,c25,c26],
[c31,c32,c33,c34,c35,c36],
[c41,c42,c43,c44,c45,c46],
[c51,c52,c53,c54,c55,c56],
[c61,c62,c63,c64,c65,c66]
];
// macro for strain and electrical field tensors
macro flucStrain(ux,uy,uz) [dx(ux),dy(uy),dz(uz),0.5*(dz(uy)+dy(uz)),0.5*(dx(uz)+dz(ux)),0.5*(dx(uy)+dy(ux))] //EOM
macro macroStrain(E11,E22,E33,E23,E13,E12)[E11,E22,E33,E23,E13,E12]//
//Define the problem
// Defines a variational problem prob where [ux, uy, uz] are the unknown displacement fields, and [vx, vy, vz] are the test functions
// int3d(mohamad): Indicates that the integration is performed over the 3D domain
// flucStrain(vx, vy, vz)': Transpose of the fluctuating strain tensor derived from the test functions.
problem prob([ux,uy,uz],[vx,vy,vz])=
int3d(Th5) ( flucStrain(vx,vy,vz)'C6x6 flucStrain(ux,uy,uz) )
+int3d(Th5)(flucStrain(vx,vy,vz)'C6x6macroStrain(E11,E22,E33,E23,E13,E12))
;
E11=1;E22=0;E33=0;E23=0;E13=0;E12=0;
prob; //resolution du probleme
// plot the variable ux
// true indicates that we want to display the numerical values of ux on the plot
// cmm stands for comment, This string will be displayed as a title or comment on the plot.
plot(ux,wait=1,value=true,cmm=“ux 1”); // affichage des valeurs de ux
plot(uy,wait=1,value=true, cmm=“uy 1”);
plot(uz,wait=1,value=true, cmm=“uz 1”);
//plot([ux,uy],wait=1,value=1,cmm=“champ de deplacement u*”); // on trace le champ de deplacement
// fluctuated values of epsilon and electric fields :
epsilonx=dx(ux);//store the computed value of the derivative of ux with respect to x, which represents the strain in the x-direction
epsilony=dy(uy); //deformation selon y
epsilonz=dz(uz);
epsilonxy=(1/2)(dx(uy)+dy(ux)); //deformation liee au cisainumbernumberllement
epsilonxz=(1/2)(dz(ux)+dy(ux));
epsilonyz=(1/2)*(dz(uy)+dy(uz));
//stresses
sigmax=c11*(epsilonx+E11)+c12*(epsilony+E22)+c13*(epsilonz+E33)+c14*(epsilonyz+E23)+c15*(epsilonxz+E13)+c16*(epsilonxy+E12);
sigmay=c21*(epsilonx+E11)+c22*(epsilony+E22)+c23*(epsilonz+E33)+c24*(epsilonyz+E23)+c25*(epsilonxz+E13)+c26*(epsilonxy+E12);
sigmaz=c31*(epsilonx+E11)+c32*(epsilony+E22)+c33*(epsilonz+E33)+c34*(epsilonyz+E23)+c35*(epsilonxz+E13)+c36*(epsilonxy+E12);
sigmayz=c41*(epsilonx+E11)+c42*(epsilony+E22)+c43*(epsilonz+E33)+c44*(epsilonyz+E23)+c45*(epsilonxz+E13)+c46*(epsilonxy+E12);
sigmaxz=c51*(epsilonx+E11)+c52*(epsilony+E22)+c53*(epsilonz+E33)+c54*(epsilonyz+E23)+c55*(epsilonxz+E13)+c56*(epsilonxy+E12);
sigmaxy=c61*(epsilonx+E11)+c62*(epsilony+E22)+c63*(epsilonz+E33)+c64*(epsilonyz+E23)+c65*(epsilonxz+E13)+c66*(epsilonxy+E12);
// define a real number chom11
// int3d(mohamad)(sigmax) computes the volume integral of sigmax over the domain mohamad. This term represents the total value of sigmax integrated over the 3D domain.
real Chom11=int3d(Th5)(sigmax)/int3d(Th5)(1.);
real Chom21=int3d(Th5)(sigmay)/int3d(Th5)(1.);
real Chom31=int3d(Th5)(sigmaz)/int3d(Th5)(1.);
real Chom41=int3d(Th5)(sigmayz)/int3d(Th5)(1.);
real Chom51=int3d(Th5)(sigmaxz)/int3d(Th5)(1.);
real Chom61=int3d(Th5)(sigmaxy)/int3d(Th5)(1.);
// Calculation of strain energy, W=0.5sigmaepsilon
wtot=0.5*((Chom11E11)+(Chom21E22)+(Chom31E33)+(Chom41E23)+(Chom51E12)+(Chom61E13));
cout<<“wtot=”<<wtot<<endl;
real Chom44 = 1-2*Chom21;
cout << “Chom11= " <<Chom11 << " Chom12= " <<Chom21 << " Chom13=” << Chom21<< “Chom14= " <<Chom41 << " Chom15= " <<Chom51 << " Chom16=” << Chom61<<endl;
cout << “Chom21= " <<Chom21 << " Chom22= " <<Chom11 << " Chom23=” << Chom21<< “Chom24= " <<Chom41 << " Chom25= " <<Chom51 << " Chom26=” << Chom61<<endl;
cout << “Chom31= " <<Chom31 << " Chom32= " <<Chom31 << " Chom33=” << Chom11<< “Chom34= " <<Chom41 << " Chom35= " <<Chom51 << " Chom36=” << Chom61<< endl;
cout << “Chom41= " <<Chom41 << " Chom42= " <<Chom41 << " Chom43=” << Chom41<< “Chom44= " <<Chom44 << " Chom45= " <<Chom51 << " Chom46=” << Chom61<<endl;
cout << “Chom51= " <<Chom51 << " Chom42= " <<Chom51 << " Chom53=” << Chom51<< “Chom54= " <<Chom41 << " Chom55= " <<Chom44 << " Chom56=” << Chom61<<endl;
cout << “Chom61= " <<Chom61 << " Chom62= " <<Chom61 << " Chom63=” << Chom61<< “Chom64= " <<Chom41 << " Chom65= " <<Chom51 << " Chom66=” << Chom44<<endl;
if(abs((test-Chom11)/test)<0.01) break;
else {
c332=c332+2129;
c112=c112-24.33333333333;
c222=c112;
c122=c122-14.33333333333;
c442=c442-4.66666666667;
c552=c442;
}
cout<<“condition1111=”<<abs((test-Chom11)/test)<<endl;
cout << "c332= " <<c332 <<endl;
cout << "c112= " <<c112 <<endl;
cout << "c222= " <<c222 <<endl;
cout << "c442= " <<c442 <<endl;
cout << "c552= " <<c552 <<endl;
cout << "R1= " <<R1 <<endl;
cout<<“wtot1111=”<<wtot<<endl;
}
if((abs(wtot-w)/w)<0.0001)break;
else
{R1=R1+0.02;
}
cout << "r1= " <<R1 <<endl;
cout << "c332= " <<c332 <<endl;
cout << "c112= " <<c112 <<endl;
cout << "c222= " <<c222 <<endl;
cout << "c442= " <<c442 <<endl;
cout << "c552= " <<c552 <<endl;
cout<<“wtot=”<<wtot<<endl;
}
cout << "R1= " <<R1 <<endl;
cout << "c332= " <<c332 <<endl;
cout << "c112= " <<c112 <<endl;
cout << "c222= " <<c222 <<endl;
cout << "c442= " <<c442 <<endl;
cout << "c552= " <<c552 <<endl;
cout<<“wtot=”<<wtot<<endl;