// // Script pour la résolution de LEWE // load "iovtk" // librairie pour la sauvegarde de la solution en format .vtu lisible avec Paraview verbosity = 0; // seulement les messages d'erreur sont affiches dans la console real tgv = 1e30; // The transducer's properties real Au = 1.; // amplitude real freq = 100e3; // frequency real Omega = 2.*pi*freq; real Sigma = 1.0/freq*0.5; real t0 = 2*Sigma; // time before emission macro f(temps) (Au*sin(Omega*(temps-t0))*exp(-(temps-t0)^2/Sigma^2)) // // Generation du maillage real hDomain = 0.03; // [m] real lDomain = 0.06; // [m] real phiPiezo = 0.0254; // [m] border B01(t = 0, lDomain){x = t; y = 0; label = 1;}; border B02(t = 0, hDomain){x = lDomain; y = t; label = 2;}; border B03(t = lDomain, (lDomain+phiPiezo)/2.){x = t; y = hDomain; label = 3;}; border B04(t = (lDomain+phiPiezo)/2., (lDomain-phiPiezo)/2.){x = t; y = hDomain; label = 4;}; // le transducteur border B05(t = (lDomain-phiPiezo)/2., 0){x = t; y = hDomain; label = 5;}; border B06(t = hDomain, 0){x = 0; y = t; label = 6;}; real lam = 2250.0/freq; int pointsPerLambda = 50; // DoF on wavelenght mesh Th = buildmesh(B01(lDomain/lam*pointsPerLambda) + B02(hDomain/lam*pointsPerLambda) + B03(((lDomain-phiPiezo)/2.)/lam*pointsPerLambda) + B04(phiPiezo/lam*pointsPerLambda) + B05(((lDomain-phiPiezo)/2.)/lam*pointsPerLambda) + B06(hDomain/lam*pointsPerLambda)); int NbTriangles = Th.nt; cout << "Nombre d'elements = " + NbTriangles << endl; // Visualisation du maillage //plot(Th, cmm="Maillage", dim=2, wait=true); //exit(1); // Generation de l'espace des elements finis fespace Vh(Th, P1); fespace Vh5(Th, [P1, P1, P1, P1, P1]); // Les parametres physiques du milieu Vh cp = 3950.0; Vh cs = 2250.0; Vh rho = 2050.0; Vh mu = rho*cs^2; Vh lambda = rho*cp^2 - 2.0*mu; real dt = 1e-8; real Tfinal = 7e-6; // temps final de propagation des ondes cout << "Pas de temps = " + dt << endl; // La forme variationelle Vh5 [Txx, Tyy, Txy, u, v]; // la solution a t+1 Vh5 [Txxc, Tyyc, Txyc, uc, vc]; // la solution a t courant Vh5 [Txxm, Tyym, Txym, um, vm]; // la solution a t-1 Vh5 [Txxt, Tyyt, Txyt, ut, vt]; // la fonction de test real tc = 0.; // le pas de temps courant Vh fc; varf vS([Txx,Tyy,Txy,u,v], [Txxt,Tyyt,Txyt,ut,vt]) = int2d(Th)([Txx,Tyy,Txy,u,v]'*[Txxt,Tyyt,Txyt,ut,vt]/dt) - int2d(Th)([-(lambda+2.0*mu)*u, -lambda*u, -mu*v, -Txx/rho, -Txy/rho]'*[dx(Txxt),dx(Tyyt),dx(Txyt),dx(ut),dx(vt)] + [-lambda*v, -(lambda+2.0*mu)*v, -mu*u, -Txy/rho, -Tyy/rho]'*[dy(Txxt),dy(Tyyt),dy(Txyt),dy(ut),dy(vt)]) + on(1,3,5, u=0,v=0) + on(2,6, u=0,v=0); varf vGamma([Txxm,Tyym,Txym,um,vm], [Txxt,Tyyt,Txyt,ut,vt]) = on(4, vm=fc); varf vM([Txxm,Tyym,Txym,um,vm], [Txxt,Tyyt,Txyt,ut,vt]) = int2d(Th)([Txxm,Tyym,Txym,um,vm]'*[Txxt,Tyyt,Txyt,ut,vt]/dt); matrix M = vM(Vh5, Vh5); matrix S = vS(Vh5, Vh5, solver=UMFPACK); // Boucle en temps int pasTemps = int(Tfinal/dt); // le nombre total de pas de temps real cpu = clock(); // le temps en seconds int[int] fforder = [1, 1, 1]; real[int] BC(Vh5.ndof); for(int i = 0; i <= pasTemps; i++){ Txxm[] = Txx[]; Txxc[] = M*Txxm[]; BC = vGamma(0, Vh5); BC = -BC; Txxc[] += BC; Txx[] = S^-1 * Txxc[]; fc = f(tc); cout << i + " " + tc + " " + (Tfinal) << endl; // Visualisation/sauvegarde solution if (!(i % 5)){ savevtk("solution_lewe_K" + NbTriangles + "_" + i + ".vtu", Th, Txx, order=fforder, dataname="Txx", bin=true); } tc += dt; } // Le temps d'execution de la boucle en temps cout << "Temps calcul = " << (clock()-cpu) << endl;