Validation test for EMT model from Inno2019

dG3D/benchmarks/Inno_EMT_test/CMakeLists.txt

+# test file
+
+set(PYFILE InnoElecMagTherm.py)
+
+set(FILES2DELETE 
+  *.msh
+  *.csv
+)
+
+add_cm3python_test(${PYFILE} "${FILES2DELETE}")

dG3D/benchmarks/Inno_EMT_test/InnoElecMagTherm.py

+#coding-Utf-8-*-
+from gmshpy import *
+
+
+#from dG3DpyDebug import*
+from dG3Dpy import*
+
+#script to launch ElecMag SMP problem with a python script
+
+import math
+
+# material law: SMP
+lawnum= 1 
+rho   = 270. # material density (Kg/m^3)
+tinitial = 273.+25. # (K)
+G=299.e6 # Shear modulus (Pa)
+Mu=0.29 # Poisson ratio
+alpha = beta=gamma=0. # parameter of anisotropy
+cp=10. # specific heat capacity (Kg m^2/K s^2) 
+Kx=Ky=Kz=237. # Thermal conductivity (W/m K)
+lx=ly=lz=1.0e4 # electrical conductivity (S/m) 
+seebeck=21.e-6
+alphaTherm=0.0 #13.e-5 # thermal dilatation
+v0=0. # initial electric potential
+E_matrix = 771.e6 # (Pa)
+mag_mu_0 = 4.e-7 * math.pi # vacuum magnetic permeability
+mag_r_smp = 20.0 # relative mag permeability of SMP
+mag_mu_x = mag_mu_y = mag_mu_z = mag_r_smp * mag_mu_0 # mag permeability smp 
+magpot0_x = magpot0_y = magpot0_z = 0.0 # initial magnetic potential
+Irms = 10. # (Ampere)
+freq = 1000. # (Hz)
+nTurnsCoil = 1000
+coilLx = coilLy = coilLz = 1. # m
+coilW = 3.e-2 # m
+
+# material law: inductor (copper coil)
+# Material parameters for copper
+lawnumind = 2
+rhoind = 8960. 
+Eind=120.e9 #youngs modulus
+Gind=48.e9 # Shear modulus
+Muind = 0.34 # Poisson ratio
+alphaind = betaind = gammaind = 0. # parameter of anisotropy
+cpind= 385.
+Kxind=Kyind=Kzind= 401. # thermal conductivity tensor components
+lxind=lyind=lzind=5.77e7 # electrical conductivity (S/m)
+seebeckind=6.5e-6
+alphaThermind= 0.0 # thermal dilatation
+v0ind = 0.
+mag_r_ind = 1.0 # relative mag permeability of inductor
+mag_mu_x_ind = mag_mu_y_ind = mag_mu_z_ind = mag_r_ind * mag_mu_0 # mag permeability inductor 
+magpot0_x_ind = magpot0_y_ind = magpot0_z_ind = 0.0 # initial magnetic potential
+Centroid_x = Centroid_y = Centroid_z = 0.0 # Centroid of inductor coil
+centralAxis_x = 0.
+centralAxis_y = 1.
+centralAxis_z = 0. # Unit vector along central axis of inductor coil
+
+# material law: vacuum/free space region
+lawnumvac = 3
+rhovac = 1.2 
+Gvac=156.e1 # Shear modulus (1.e-5 less than SMP)
+Muvac = 0.35 # Poisson ratio (same as SMP)
+Evac= 2.0 * Gvac * (1.0 + Muvac) #youngs modulus
+alphavac = betavac = gammavac = 0. # parameter of anisotropy
+cpvac= 1012.
+Kxvac=Kyvac=Kzvac= 26.e-3 # thermal conductivity tensor components
+lxvac=lyvac=lzvac= 1.e-12 # electrical conductivity
+seebeckvac= 0.
+alphaThermvac= 0.0 # thermal dilatation
+v0vac = 0.
+mag_r_vac = 1.0 # relative mag permeability of inductor
+mag_mu_x_vac = mag_mu_y_vac = mag_mu_z_vac = mag_r_vac * mag_mu_0 # mag permeability inductor 
+magpot0_x_vac = magpot0_y_vac = magpot0_z_vac = 0.0 # initial magnetic potential
+
+# geometry
+geofile="smp.geo" # name of mesh file
+meshfile="smp.msh" # name of mesh file
+
+# solver
+sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
+soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
+nstep = 50 # number of step (used only if soltype=1)
+ftime =1.0/freq # Final time (used only if soltype=1)
+tol=1.e-3   # relative tolerance for NR scheme (used only if soltype=1)
+nstepArch=1 # Number of step between 2 archiving (used only if soltype=1)
+fullDg = 0 #O = CG, 1 = DG
+beta1 = 1000.
+eqRatio =1.e8
+
+#  compute solution and BC (given directly to the solver
+# creation of law
+
+# material law for smp
+lawsmp = LinearElecMagTherMechDG3DMaterialLaw(lawnum,rho,E_matrix,E_matrix,E_matrix,Mu,Mu,Mu,G,G,G,alpha,beta,gamma,
+tinitial,Kx,Ky,Kz,alphaTherm,alphaTherm,alphaTherm,lx,ly,lz,seebeck,cp,v0, mag_mu_x, mag_mu_y, mag_mu_z, magpot0_x,
+magpot0_y, magpot0_z, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW)
+
+#lawsmp.setUseBarF(True)
+
+# material law for inductor region
+lawind = LinearElecMagInductorDG3DMaterialLaw(lawnumind,rhoind,Eind,Eind,Eind,Muind,Muind,Muind,Gind,Gind,Gind,alphaind,betaind,gammaind,
+tinitial,Kxind,Kyind,Kzind,alphaThermind,alphaThermind,alphaThermind,lxind,lyind,lzind,seebeckind,cpind,v0ind, mag_mu_x_ind, mag_mu_y_ind, mag_mu_z_ind, 
+magpot0_x_ind,magpot0_y_ind, magpot0_z_ind, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW, 
+Centroid_x,Centroid_y,Centroid_z,centralAxis_x,centralAxis_y,centralAxis_z)
+
+#lawind.setUseBarF(True)
+
+#material law for free space (vacuum)
+lawvac = LinearElecMagTherMechDG3DMaterialLaw(lawnumvac,rhovac,Evac,Evac,Evac,Muvac,Muvac,Muvac,Gvac,Gvac,Gvac,alphavac,betavac,gammavac,tinitial,
+Kxvac,Kyvac,Kzvac,alphaThermvac,alphaThermvac,alphaThermvac,lxvac,lyvac,lzvac,seebeckvac,cpvac,v0vac, mag_mu_x_vac, mag_mu_y_vac, mag_mu_z_vac, magpot0_x_vac,
+magpot0_y_vac, magpot0_z_vac, Irms, freq, nTurnsCoil, coilLx, coilLy, coilLz, coilW)
+
+#lawvac.setUseBarF(True)
+
+# creation of ElasticField
+SMPfield = 1000 # number of the field (physical number of SMP)
+Indfield = 2000 # number of the field (physical number of Inductor region)
+Vacfield = 3000 # number of the field (physical number of Vacuum)
+SurfVacOut = 3333 # number of the field (outer surface of Vacuum) Used for BC
+SurfInd = 2222 # outer surface of inductor region
+SurfSmp = 11110 # outer surface of conductor SMP
+SmpRef=11111 # top face of outer surface of SMP
+SmpRefBottom=11115 # bottom face of outer surface of SMP
+SmpRefPoint=11112 # single point on top face of SMP
+
+space1 = 0 # function space (Lagrange=0)
+
+SMP_field = ElecMagTherMechDG3DDomain(10,SMPfield,space1,lawnum,fullDg,eqRatio,3)
+Inductor_field = ElecMagTherMechDG3DDomain(10,Indfield,space1,lawnumind,fullDg,eqRatio,3)
+Vacuum_field = ElecMagTherMechDG3DDomain(10,Vacfield,space1,lawnumvac,fullDg,eqRatio,3)
+
+SMP_field.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
+Inductor_field.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
+Vacuum_field.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
+SMP_field.setConstitutiveCurlEqRatio(eqRatio)
+Inductor_field.setConstitutiveCurlEqRatio(eqRatio)
+Vacuum_field.setConstitutiveCurlEqRatio(eqRatio)
+SMP_field.stabilityParameters(beta1)
+Inductor_field.stabilityParameters(beta1)
+Vacuum_field.stabilityParameters(beta1)
+#Matixfield.setConstitutiveExtraDofDiffusionStabilityParameters(beta1,bool(1))
+#Fiberfield.setConstitutiveExtraDofDiffusionStabilityParameters(beta1,bool(1))
+
+thermalSource=True
+mecaSource=False
+
+SMP_field.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+Vacuum_field.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+Inductor_field.setConstitutiveExtraDofDiffusionAccountSource(thermalSource,mecaSource)
+
+#interface // we need to create the ininterdomain after adding the domains
+#myinterfield=interDomainBetween3D(12,Matixfield,Fiberfield) 
+#myinterfield.setConstitutiveExtraDofDiffusionEqRatio(eqRatio)
+#myinterfield.stabilityParameters(beta1)
+#myinterfield.setConstitutiveExtraDofDiffusionStabilityParameters(beta1,bool(1))
+#myinterfield.matrixByPerturbation(0,1e-8)
+
+#SMP_field.matrixByPerturbation(1,1,1,1.e-3)
+#Inductor_field.matrixByPerturbation(1,1,1,1.e-3)
+#Vacuum_field.matrixByPerturbation(1,1,1,1.e-3)
+
+# creation of Solver
+mysolver = nonLinearMechSolver(1000)
+#mysolver.createModel(geofile,meshfile,3,1)
+mysolver.loadModel(meshfile)
+mysolver.addDomain(SMP_field)
+mysolver.addDomain(Inductor_field)
+mysolver.addDomain(Vacuum_field)
+mysolver.addMaterialLaw(lawsmp)
+mysolver.addMaterialLaw(lawind)
+mysolver.addMaterialLaw(lawvac)
+mysolver.Scheme(soltype)
+mysolver.Solver(sol)
+mysolver.snlData(nstep,ftime,tol,1.e-9)#forsmp
+mysolver.options("-ksp_type preonly")
+mysolver.options("-pc_type lu")
+
+#mysolver.snlData(nstep,ftime,tol,1.e-19)
+mysolver.stepBetweenArchiving(nstepArch)
+mysolver.snlManageTimeStep(25, 3, 2, 10)
+
+# BC
+##cyclicFunctionDisp=cycleFunctionTime(0., 0., ftime/5., -0.0001, ftime , -0.0001);
+#mysolver.displacementBC("Face",83,2,0.)
+mysolver.displacementBC("Volume",SMPfield,0,0.0)
+mysolver.displacementBC("Volume",SMPfield,1,0.0)
+mysolver.displacementBC("Volume",SMPfield,2,0.0)
+mysolver.displacementBC("Volume",Indfield,0,0.0)
+mysolver.displacementBC("Volume",Indfield,1,0.0)
+mysolver.displacementBC("Volume",Indfield,2,0.0)
+mysolver.displacementBC("Volume",Vacfield,0,0.0)
+mysolver.displacementBC("Volume",Vacfield,1,0.0)
+mysolver.displacementBC("Volume",Vacfield,2,0.0)
+
+#thermal BC
+cyclicFunctionTemp1=cycleFunctionTime(0., tinitial,ftime,tinitial);
+#mysolver.displacementBC("Face",85,3,cyclicFunctionTemp1)
+#mysolver.initialBC("Volume","Position",Mfield,3,tinitial)
+mysolver.initialBC("Volume","Position",SMPfield,3,tinitial)
+mysolver.initialBC("Volume","Position",Indfield,3,tinitial)
+mysolver.initialBC("Volume","Position",Vacfield,3,tinitial)
+mysolver.displacementBC("Volume",Indfield,3,cyclicFunctionTemp1)
+#mysolver.displacementBC("Volume",SMPfield,3,cyclicFunctionTemp1)
+#mysolver.displacementBC("Volume",Vacfield,3,cyclicFunctionTemp1)
+mysolver.displacementBC("Face",SurfVacOut,3,cyclicFunctionTemp1)
+#mysolver.displacementBC("Face",SmpRef,3,cyclicFunctionTemp1)
+#mysolver.displacementBC("Face",SmpRefBottom,3,cyclicFunctionTemp1)
+
+#electrical BC
+#cyclicFunctionvolt1=cycleFunctionTime(0., 0.,  ftime,10.);
+#mysolver.displacementBC("Face",85,4,cyclicFunctionvolt1)
+mysolver.displacementBC("Volume",Indfield,4,0.0)
+mysolver.initialBC("Volume","Position",SMPfield,4,0.0)
+mysolver.initialBC("Volume","Position",Vacfield,4,0.0)
+mysolver.displacementBC("Face",SurfVacOut,4,0.0)
+#mysolver.displacementBC("Face",SmpRef,4,0.0)
+mysolver.displacementBC("Node",SmpRefPoint,4,0.0)
+
+#magentic
+mysolver.curlDisplacementBC("Face",SurfVacOut,5,0.0) # comp may also be 5
+#Gauging the edges using tree-cotree methodcd
+PhysicalCurves = "" 		# input required as a string of comma separated ints
+PhysicalSurfaces = "11110,3333" # input required as a string of comma separated ints
+PhysicalVolumes = "1000,2000,3000" 	# input required as a string of comma separated ints
+OutputPhysical = 55 		# input required as a int
+mysolver.createTreeForBC(PhysicalCurves,PhysicalSurfaces,PhysicalVolumes,OutputPhysical)
+mysolver.curlDisplacementBC("Edge",OutputPhysical,5,0.0)
+
+#mysolver.internalPointBuildView("svm",IPField.SVM, 1, 1)
+#mysolver.internalPointBuildView("sig_xx",IPField.SIG_XX, 1, 1)
+#mysolver.internalPointBuildView("sig_yy",IPField.SIG_YY, 1, 1)
+#mysolver.internalPointBuildView("sig_zz",IPField.SIG_ZZ, 1, 1)
+#mysolver.internalPointBuildView("sig_xy",IPField.SIG_XY, 1, 1)
+#mysolver.internalPointBuildView("sig_yz",IPField.SIG_YZ, 1, 1)
+#mysolver.internalPointBuildView("sig_xz",IPField.SIG_XZ, 1, 1)
+mysolver.internalPointBuildView("temperature",IPField.TEMPERATURE, 1, 1)
+mysolver.internalPointBuildView("jyx",IPField.THERMALFLUX_X, 1, 1)
+mysolver.internalPointBuildView("jyy",IPField.THERMALFLUX_Y, 1, 1)
+mysolver.internalPointBuildView("jyz",IPField.THERMALFLUX_Z, 1, 1)
+mysolver.internalPointBuildView("voltage",IPField.VOLTAGE, 1, 1)
+mysolver.internalPointBuildView("jex",IPField.ELECTRICALFLUX_X, 1, 1)
+mysolver.internalPointBuildView("jey",IPField.ELECTRICALFLUX_Y, 1, 1)
+mysolver.internalPointBuildView("jez",IPField.ELECTRICALFLUX_Z, 1, 1)
+#mysolver.internalPointBuildView("ax",IPField.MAGNETICVECTORPOTENTIAL_X, 1, 1)
+#mysolver.internalPointBuildView("ay",IPField.MAGNETICVECTORPOTENTIAL_Y, 1, 1)
+#mysolver.internalPointBuildView("az",IPField.MAGNETICVECTORPOTENTIAL_Z, 1, 1)
+mysolver.internalPointBuildView("bx",IPField.MAGNETICINDUCTION_X, 1, 1)
+mysolver.internalPointBuildView("by",IPField.MAGNETICINDUCTION_Y, 1, 1)
+mysolver.internalPointBuildView("bz",IPField.MAGNETICINDUCTION_Z, 1, 1)
+mysolver.internalPointBuildView("hx",IPField.MAGNETICFIELD_X, 1, 1)
+mysolver.internalPointBuildView("hy",IPField.MAGNETICFIELD_Y, 1, 1)
+mysolver.internalPointBuildView("hz",IPField.MAGNETICFIELD_Z, 1, 1)
+mysolver.internalPointBuildView("js0_x",IPField.INDUCTORSOURCEVECTORFIELD_X, 1, 1)
+mysolver.internalPointBuildView("js0_y",IPField.INDUCTORSOURCEVECTORFIELD_Y, 1, 1)
+mysolver.internalPointBuildView("js0_z",IPField.INDUCTORSOURCEVECTORFIELD_Z, 1, 1)
+mysolver.internalPointBuildView("wAV_x",IPField.EMSOURCEVECTORFIELD_X, 1, 1)
+mysolver.internalPointBuildView("wAV_y",IPField.EMSOURCEVECTORFIELD_Y, 1, 1)
+mysolver.internalPointBuildView("wAV_z",IPField.EMSOURCEVECTORFIELD_Z, 1, 1)
+
+#mysolver.archivingIPOnPhysicalGroup("Volume",Mfield, IPField.SIG_ZZ,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Mfield, IPField.STRAIN_ZZ,IPField.MEAN_VALUE,nstepArch);
+mysolver.archivingIPOnPhysicalGroup("Volume",SMPfield, IPField.TEMPERATURE,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Mfield, IPField.VOLTAGE,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Ffield, IPField.SIG_ZZ,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Ffield, IPField.STRAIN_ZZ,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Ffield, IPField.TEMPERATURE,IPField.MEAN_VALUE,nstepArch);
+mysolver.archivingIPOnPhysicalGroup("Volume",SMPfield, IPField.VOLTAGE,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Indfield, IPField.VOLTAGE,IPField.MEAN_VALUE,nstepArch);
+#mysolver.archivingIPOnPhysicalGroup("Volume",Vacfield, IPField.VOLTAGE,IPField.MEAN_VALUE,nstepArch);
+
+mysolver.archivingForceOnPhysicalGroup("Face", SurfSmp, 3, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfSmp, 4, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfInd, 4, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfVacOut, 4, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfSmp, 5, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfInd, 5, 1);
+mysolver.archivingForceOnPhysicalGroup("Face", SurfVacOut, 5, 1);
+
+mysolver.archivingNodeDisplacement(SmpRefPoint,0,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,1,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,2,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,3,1);
+mysolver.archivingNodeDisplacement(SmpRefPoint,4,1);
+
+#mysolver.archivingNodeIP(SmpRef, IPField.THERMALFLUX_X,IPField.MEAN_VALUE,1);
+#mysolver.archivingNodeIP(SmpRef, IPField.ELECTRICALFLUX_X,IPField.MEAN_VALUE,1);
+
+mysolver.solve()
+
+check = TestCheck()
+check.equal(-6.995019e+11,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSmp, 3),1.e-5)
+check.equal(0.000000e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSmp, 4),1.e-5)
+check.equal(6.181238e-05,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfInd, 4),1.e-5)
+check.equal(1.826831e-08,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfVacOut, 4),1.e-5)
+check.equal(2.529895e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfSmp, 5),1.e-5)
+check.equal(-2.498196e+00,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfInd, 5),1.e-5)
+check.equal(-3.522163e+02,mysolver.getArchivedForceOnPhysicalGroup("Face", SurfVacOut, 5),1.e-5)
+check.equal(0.0,mysolver.getArchivedNodalValue(SmpRefPoint,0,mysolver.displacement),1.e-5)
+check.equal(0.0,mysolver.getArchivedNodalValue(SmpRefPoint,1,mysolver.displacement),1.e-5)
+check.equal(0.0,mysolver.getArchivedNodalValue(SmpRefPoint,2,mysolver.displacement),1.e-5)
+check.equal(2.980000e+02,mysolver.getArchivedNodalValue(SmpRefPoint,3,mysolver.displacement),1.e-5)
+check.equal(-3.511321e-04,mysolver.getArchivedNodalValue(SmpRefPoint,4,mysolver.displacement),1.e-5)
+

dG3D/benchmarks/Inno_EMT_test/hollow_smp.geo

+Include "hollow_smp_data.geo";
+
+SetFactory("OpenCASCADE");
+
+Mesh.Optimize = 1;
+
+// circular x-section coil
+vC()+=newv; Cylinder(newv) = {0,-Ly/2,0, 0,Ly,0, RintCoil};
+vC()+=newv; Cylinder(newv) = {0,-Ly/2,0, 0,Ly,0, RextCoil};
+
+vCoil()+=newv; BooleanDifference(newv) = { Volume{vC(1)}; Delete; }{ Volume{vC(0)}; Delete; };
+
+// Smp 
+//vCoreE()=newv; Cylinder(newv) = { 0,-H/2,0,  0,H,0,  R};
+vE()+=newv; Cylinder(newv) = {0,-H/2,0, 0,H,0, RintSMP};
+vE()+=newv; Cylinder(newv) = {0,-H/2,0, 0,H,0, RextSMP};
+
+vCoreE()+=newv; BooleanDifference(newv) = { Volume{vE(1)}; Delete; }{ Volume{vE(0)}; Delete; };
+
+// Air around
+vA()+=newv; Sphere(newv) = {0,0,0, Rint};
+vA()+=newv; Sphere(newv) = {0,0,0, Rext};
+
+vAir()+=newv; BooleanDifference(newv) = { Volume{vA(1)}; Delete; }{ Volume{vA(0)}; };
+vAir()+=newv; BooleanDifference(newv) = { Volume{vA(0)}; Delete; }{ Volume{vCoreE(),vCoil()}; };
+
+If(Flag_Symmetry)
+  xa =  -Rext;   ya =  -Rext;  za =  -Rext;
+  dxa =  2*Rext; dya = 2*Rext; dza = 2*Rext;
+
+  If(Flag_Symmetry==1)
+    vAux=newv; Box(newv) = {xa, ya, 0*za, dxa, dya, -dza/2};
+  EndIf
+  If(Flag_Symmetry==2)
+    vAux=newv; Box(newv) = {0*xa, ya, 0*za, dxa/2, dya, -dza/2};
+  EndIf
+
+  For k In {0:#vCoil()-1}
+    vvv=newv; BooleanIntersection(newv) = { Volume{vCoil(k)}; Delete; }{ Volume{vAux}; };
+    vCoil(k) = vvv;
+  EndFor
+  For k In {0:#vAir()-1}
+    vvv=newv; BooleanIntersection(newv) = { Volume{vAir(k)}; Delete; }{ Volume{vAux}; };
+    vAir(k) = vvv;
+  EndFor
+
+  For k In {0:#vCoreE()-1}
+    vvv=newv; BooleanIntersection(newv) = { Volume{vCoreE(k)}; Delete; }{ Volume{vAux}; };
+    vCoreE(k) = vvv;
+  EndFor
+
+  //vvv=newv; BooleanIntersection(newv) = { Volume{vCoreE(0)}; Delete; }{ Volume{vAux}; Delete; };
+  //vCoreE(0) = vvv;
+EndIf
+
+BooleanFragments{ Volume{vAir(), vCoreE(), vCoil()}; Delete; }{} // This needs to be done at the end
+
+// Adapting mesh size
+// characteristic lengths
+lc0 = Pi*Rint/5; // air
+lc1  = wcoil/2; // coil 
+lc2  = RextSMP/6; // smp
+
+Characteristic Length { PointsOf{ Volume{vAir(0)}; } } = lc0;
+Characteristic Length { PointsOf{ Volume{vCoil()}; } } = lc1;
+Characteristic Length { PointsOf{ Volume{vCoreE()}; } } = lc2;
+
+
+// Boundary conditions
+tol = 1e-5;
+cut_xy() = Surface In BoundingBox {-Rext-tol,-Rext-tol,-tol, 2*(Rext+tol), 2*(Rext+tol), 2*tol}; // 1/2, 1/4
+cut_yz() = Surface In BoundingBox {-tol,-Rext-tol,-Rext-tol, 2*tol, 2*(Rext+tol), 2*(Rext+tol)}; // 1/4
+
+all_vol() = Volume '*';
+Printf("all_vol()=", all_vol());
+bndAir() = CombinedBoundary{Volume{all_vol()};};
+bndAir() -= {cut_xy(),cut_yz()};
+
+
+//=================================================
+// Some colors... for aesthetics :-)
+//=================================================
+Recursive Color SkyBlue {Volume{vAir()};}
+Recursive Color SteelBlue {Volume{vCoreE()};}
+Recursive Color Red {Volume{vCoil()};}
+
+//=================================================
+// Physical regions for FE analysis with GetDP
+//=================================================
+
+Physical Volume(ECORE) = vCoreE();
+bnd_vCoreE() = CombinedBoundary{Volume{vCoreE()};};
+bnd_vCoreE() -= {cut_xy(),cut_yz()};
+Physical Surface(SKINECORE) = bnd_vCoreE();
+
+Physical Surface(REFCORE)= bnd_vCoreE(1);
+Physical Surface(REFCOREBOTTOM)= bnd_vCoreE(2);
+pts() = PointsOf {Surface { bnd_vCoreE(1) }; };
+Physical Point(REFCOREPOINT)= pts(0);
+
+Physical Volume(COIL) = vCoil();
+bnd_vCoil() = CombinedBoundary{Volume{vCoil()};};
+bnd_vCoil() -= {cut_xy(),cut_yz()};
+Physical Surface(SKINCOIL) = bnd_vCoil();
+
+If(Flag_Infinity==0)
+  Physical Volume(AIR) = {vAir()};
+EndIf
+If(Flag_Infinity==1)
+  Physical Volume(AIR) = vAir(1);
+  Physical Volume(AIRINF) = vAir(0);
+EndIf
+Physical Surface(SURF_AIROUT) = bndAir();
+
+Physical Surface(CUT_XY) = {cut_xy()};
+Physical Surface(CUT_YZ) = {cut_yz()}; // BC if symmetry

dG3D/benchmarks/Inno_EMT_test/hollow_smp_data.geo

+// Geometrical data for inductor model
+
+cm = 1e-2; // Unit
+
+pp  = "Input/10Geometric dimensions/0";
+pp2 = "Input/10Geometric dimensions/1Shell radius/";
+ppm = "Input/11Mesh control (Nbr of divisions)/";
+
+DefineConstant[
+   Flag_3Dmodel = 1
+   Flag_boolean = 1
+   Flag_Symmetry = {0, Choices{0="Full",1="Half",2="One fourth"},
+					Name "Input/01Symmetry type", Highlight "Blue"},
+   Flag_Infinity = {0, Choices{0,1},
+					Name "Input/01Use shell transformation to infinity"}
+];
+
+close_menu = 1;
+colorro  = "LightGrey";
+colorpp = "Ivory";
+
+
+// Smp
+RintSMP = 0.5*cm;
+RextSMP = 1*cm;
+H = 2*cm;
+
+//Coil
+RintCoil = 4.5*cm;
+RextCoil = 7.5*cm;
+
+DefineConstant[
+  //Lx    = {9*cm, Name StrCat[pp, "0Coil length along x-axis [m]"], Highlight Str[colorpp]}
+  Ly    = {100*cm, Name StrCat[pp, "0Coil length along y-axis [m]"], Highlight Str[colorpp]}
+  //Lz    = {9*cm, Name StrCat[pp, "2Coil length along z-axis [m]"], Highlight Str[colorpp]}
+  wcoil = {3*cm, Name StrCat[pp, "1Coil width [m]"], Highlight Str[colorpp]}
+  RintCoil = {4.5*cm, Name StrCat[pp, "2Coil inner radius [m]"], Highlight Str[colorpp]}
+  RextCoil = {7.5*cm, Name StrCat[pp, "3Coil outer radius [m]"], Highlight Str[colorpp]}
+  RintSMP  = {0.5*cm, Name StrCat[pp, "4Sample inner radius [m]"], Highlight Str[colorpp]}
+  RextSMP  = {1*cm, Name StrCat[pp, "5Sample outer radius [m]"], Highlight Str[colorpp]}
+  H     = {2*cm, Name StrCat[pp, "6Sample height [m]"], Highlight Str[colorpp]}
+];
+
+// radious for surrounding air with transformation to infinity
+If(Flag_Infinity==1)
+  label_Rext = "1Outer [m]";
+Else
+  label_Rext = "1[m]";
+EndIf
+
+DefineConstant[
+  Rint = {200*cm, Min 0.15, Max 0.9, Step 0.1, Name StrCat[pp2, "0Inner [m]"],
+    Visible (Flag_Infinity==1), Highlight Str[colorpp] },
+  Rext = {280*cm, Min Rint, Max 1, Step 0.1, Name StrCat[pp2, StrCat["1", label_Rext]],
+    Visible 1, Highlight Str[colorpp] }
+];
+
+Val_Rint = Rint;
+Val_Rext = Rext;
+
+
+IA = 10;
+Nw = 1000;
+
+sigma_al = 3.72e7 ; // conductivity of aluminum [S/m]
+sigma_cu = 5.77e7  ; // conductivity of copper [S/m]
+
+// ----------------------------------------------------
+// Numbers for physical regions in .geo and .pro files
+// ----------------------------------------------------
+
+ECORE = 1000;
+SKINECORE = 11110;
+REFCORE = 11111;
+REFCOREBOTTOM = 11115;
+REFCOREPOINT = 11112;
+
+COIL  = 2000;
+SKINCOIL = 2222;
+SKINCOIL_ = 2223;
+
+SURF_ELEC0 = 2333;
+SURF_ELEC1 = 2444;
+CUTCOIL = 2555;
+
+AIR    = 3000;
+AIRINF = 3100;
+AIRGAP = 3200;
+
+// Lines and surfaces for boundary conditions
+SURF_AIROUT = 3333;
+
+AXIS_Y = 10000;
+CUT_YZ = 11000;
+CUT_XY = 12000;