load "msh3" load "tetgen" load "mshmet" load "medit" load "iovtk"; load "MUMPS_seq" load "gmsh" // IOVTK ARRAY: per salvare i risultati in formato vtk e visualizzarli successivamente con Paraview int[int] ffordervel=[1]; // Leggere la mesh iniziale // Se letta da file: //mesh3 Th3("./mesh/initial.mesh"); //meshS ThS=gmshloadS("CoreFull-NoGrid.msh"); //plot(ThS, cmm="ThS"); //VA IN CRASH SUBITO // Se creata direttamente nel main, ad esempio //int nn = 30; //int[int] labs = [1, 2, 2, 1, 1, 2]; // Label numbering //mesh3 Th3 = cube(nn, nn, nn, label=labs); real cx1 = 0.0; // cerchio esterno real cy1 = 0.0; real r1 = 16.9; real cx2 = 0.0; // canale centrale real cy2 = 0.0; real r2 = 2.5; real cx3 = 8.0; // reg real cy3 = 5.0; real r3 = 2.5; real cx4 = 0.0; // shim real cy4 = -8.0; real r4 = 2.5; real cx5 = -10.0; // trans real cy5 = +8.0; real r5 = 2.5; real zmin=0,zmax=12.9; int[int] rup=[0,1], rdown=[0,2], rmid=[1,3]; border C1(t1=0,2*pi){x=cx1+r1*cos(t1);y=cy1+r1*sin(t1);label=1;}; // cerchio esterno border C2(t1=0,2*pi){x=cx2+r2*cos(t1);y=cy2+r2*sin(t1);label=2;}; // canale centrale border C3(t1=0,2*pi){x=cx3+r3*cos(t1);y=cy3+r3*sin(t1);label=3;}; // reg border C4(t1=0,2*pi){x=cx4+r4*cos(t1);y=cy4+r4*sin(t1);label=4;}; // shim border C5(t1=0,2*pi){x=cx5+r5*cos(t1);y=cy5+r5*sin(t1);label=5;}; // trans int m=5; mesh Thcercle = buildmesh(C1(20*m)+C2(+5*m)+C3(+5*m)+C4(+5*m)+C5(+5*m)); plot(Thcercle,wait =1); mesh3 Th3=buildlayers(Thcercle,coef=2, 4*m, zbound=[zmin,zmax], labelmid=rmid, labelup = rup, labeldown = rdown); // BOUNDARY CONDITIONS // ... // FUNCTION SPACE fespace Xh0(Th3, P0); fespace Xh1(Th3, P1); // METRIC Xh1 m11, m12, m22, m13, m23, m33; // ELEMENT VOLUME Xh0 vol; varf vvol(unused, chiK) = int3d(Th3)(chiK); vol[] = vvol(0, Xh0); // ADAPTATION TOLERANCE real tau = 2; // ADAPTATION CONSTANTS real hmax = 1.; real hmin = 1e-2; int nadapt = 10; real toll = 15/100., err = 150+toll, ne; include "makemetrica.edp"; ne = Th3.nt; Xh1 v = x*y*z; savevtk("nome_file_input.vtk", Th3, v, dataname="Function", order=ffordervel); // ADAPTATION cout << "Adapting..." << endl; makemetrica(tau, v[], hmax, hmin); savemesh(Th3, "mymesh.mesh"); savesol("mysol.sol", Th3, [m11, m12, m22, m13, m23, m33]); exec("mmg3d mymesh.mesh -sol mysol.sol -m 10000"); mesh3 Th3b("mymesh.o.mesh"); cout << "Number of elements = " << Th3b.nt << endl; Th3 = Th3b; // INTERPOLATION OF ALL THE VARIABLES ON THE NEW GRID // ... v = v; err = abs(Th3.nt - ne)/ne; cout << "err = " << err << endl; savevtk("nome_file_output.vtk", Th3, v, dataname="Function", order=ffordervel);