More generalized way of generating new merged mesh for weak coupling test,...

More generalized way of generating new merged mesh for weak coupling test, change of eqRatio in strong coupling test

dG3D/benchmarks/linearElectroMagnetic/LinearElecMag.py

@@ -85,7 +85,7 @@ sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
 soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
 nstep = 30 # 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)
+tol=1.e-6   # 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.
@@ -223,7 +223,7 @@ 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)
+mysolver.displacementBC("Node",SmpRefPoint,4,0.0)
 
 #magentic
 mysolver.curlDisplacementBC("Face",SurfVacOut,5,0.0) # comp may also be 5

dG3D/benchmarks/strongCouplingSMP/StrongCouplingEMTM.py

@@ -82,14 +82,14 @@ meshfile="undeformedSMPFullDomain.msh" # name of mesh file
 sol = 2  # Gmm=0 (default) Taucs=1 PETsc=2
 soltype = 1 # StaticLinear=0 (default) StaticNonLinear=1
 #EMncycles = 5; # num of sinusoidal cycles of EM frequency application
-nstep = 80 # 8*EMncycles # number of step (used only if soltype=1)
+nstep = 200 # 8*EMncycles # number of step (used only if soltype=1)
 ftime = 5.e-3 # EMncycles/freq; # according to characteristic time
 #EMftime = EMncycles/freq # Final time (used only if soltype=1)
 tol=1.e-6   # 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
+eqRatio =1.e0
 
 #  compute solution and BC (given directly to the solver
 # creation of law

dG3D/benchmarks/weakCouplingSMP/TM_weakCoupling.py

@@ -149,7 +149,7 @@ mysolver.displacementBC("Volume",SMPfield,4,0.0)
 mysolver.curlDisplacementBC("Volume",SMPfield,5,0.0) # comp may also be 5
 #Gauging the edges using tree-cotree method
 PhysicalCurves = "" 		# input required as a string of comma separated ints
-PhysicalSurfaces = "11110,11111,11112,11115" # input required as a string of comma separated ints
+PhysicalSurfaces = "11110" # input required as a string of comma separated ints
 PhysicalVolumes = "1000" 	# input required as a string of comma separated ints
 OutputPhysical = 55 		# input required as a int
 mysolver.createTreeForBC(PhysicalCurves,PhysicalSurfaces,PhysicalVolumes,OutputPhysical)

dG3D/benchmarks/weakCouplingSMP/new_merge_mesh.geo

+//Merge "deformed_mesh_2.msh";
+Include "deformed_mesh_filename.geo";
+sl1 = newsl;
+Surface Loop(sl1) = Surface{:};
+//+
+
+cm = 1e-2; // Unit
+pp  = "Input/10Geometric dimensions/0";
+pp2 = "Input/10Geometric dimensions/1Shell radius/";
+colorpp = "Ivory";
+
+//Coil
+//+
+SetFactory("Built-in");
+//+
+p1 = newp;
+Point(p1) = {0, 0.5, 0, 1.0};
+//+
+p2 = newp;
+Point(p2) = {0.045, 0.5, 0, 1.0};
+//+
+p3 = newp;
+Point(p3) = {-0.045, 0.5, 0, 1.0};
+//+
+p4 = newp;
+Point(p4) = {0, 0.5, 0.045, 1.0};
+//+
+p5 = newp;
+Point(p5) = {0, 0.5, -0.045, 1.0};
+//+
+c1 = newc;
+Circle(c1) = {p5, p1, p3};
+//+
+c2 = newc;
+Circle(c2) = {p3, p1, p4};
+//+
+c3 = newc;
+Circle(c3) = {p4, p1, p2};
+//+
+c4 = newc;
+Circle(c4) = {p2, p1, p5};
+//+
+p6 = newp;
+Point(p6) = {0.075, 0.5, 0, 1.0};
+//+
+p7 = newp;
+Point(p7) = {-0.075, 0.5, 0, 1.0};
+//+
+p8 = newp;
+Point(p8) = {0, 0.5, 0.075, 1.0};
+//+
+p9 = newp;
+Point(p9) = {0, 0.5, -0.075, 1.0};
+//+
+c5 = newc;
+Circle(c5) = {p9, p1, p7};
+//+
+c6 = newc;
+Circle(c6) = {p7, p1, p8};
+//+
+c7 = newc;
+Circle(c7) = {p8, p1, p6};
+//+
+c8 = newc;
+Circle(c8) = {p6, p1, p9};
+//+
+cl1 = newcl;
+Curve Loop(cl1) = {c5,c6,c7,c8};
+//+
+cl2 = newcl;
+Curve Loop(cl2) = {c1,c2,c3,c4};
+//+
+s1 = news;
+Surface(s1) = {cl1, cl2};
+//+
+
+Transfinite Curve {c1,c2,c3,c4} = 11 Using Progression 1;
+//+
+Transfinite Curve {c5,c6,c7,c8} = 11 Using Progression 1;
+
+//+
+e() = Extrude {0, -1, 0} {
+  Surface{s1}; Layers {10};
+};
+
+/*
+Printf("top curve = %g", e(0));
+Printf("vol = %g", e(1));
+Printf("surf = %g", e(2));
+Printf("surf = %g", e(3));
+Printf("surf = %g", e(4));
+Printf("surf = %g", e(5));
+Printf("surf = %g", e(6));
+Printf("surf = %g", e(7));
+Printf("surf = %g", e(8));
+Printf("surf = %g", e(9));
+*/
+
+Transfinite Surface {e(2),e(3),e(4),e(5),e(6),e(7),e(8),e(9)};
+
+sl2 = newsl;
+Surface Loop(sl2) = Boundary{Volume{e(1)}; };
+
+//+
+Physical Surface(2222) = Boundary{Volume{e(1)}; };
+Physical Volume(2000) = {e(1)};
+
+DefineConstant[
+  Ly    = {100*cm, Name StrCat[pp, "3Coil length along y-axis [m]"], Highlight Str[colorpp]}
+  wcoil = {3*cm, Name StrCat[pp, "4Coil width [m]"], Highlight Str[colorpp]}
+  RintCoil = {4.5*cm, Name StrCat[pp, "5Coil inner radius [m]"], Highlight Str[colorpp]}
+  RextCoil = {7.5*cm, Name StrCat[pp, "6Coil outer radius [m]"], Highlight Str[colorpp]}
+];
+
+//+
+p10 = newp;
+Point(p10) = {0, 0, 0, 1.0};
+//+
+
+//Air
+//+
+p11 = newp;
+Point(p11) = {2.8, 0, 0, 1.0};
+//+
+p12 = newp;
+Point(p12) = {-2.8, 0, 0, 1.0};
+//+
+p13 = newp;
+Point(p13) = {0, 2.8, 0, 1.0};
+//+
+p14 = newp;
+Point(p14) = {0, -2.8, 0, 1.0};
+//+
+c9 = newc;
+Circle(c9) = {p11, p10, p13};
+//+
+c10 = newc;
+Circle(c10) = {p13, p10, p12};
+//+
+c11 = newc;
+Circle(c11) = {p12, p10, p14};
+//+
+c12 = newc;
+Circle(c12) = {p14, p10, p11};
+//+
+
+
+//+
+p15 = newp;
+Point(p15) = {0, 0, 2.8, 1.0};
+//+
+p16 = newp;
+Point(p16) = {0, 0, -2.8, 1.0};
+//+
+c13 = newc;
+Circle(c13) = {p12, p10, p15};
+//+
+c14 = newc;
+Circle(c14) = {p15, p10, p11};
+//+
+c15 = newc;
+Circle(c15) = {p11, p10, p16};
+//+
+c16 = newc;
+Circle(c16) = {p16, p10, p12};
+//+
+
+
+
+//+
+c17 = newc;
+Circle(c17) = {p13, p10, p15};
+//+
+c18 = newc;
+Circle(c18) = {p15, p10, p14};
+//+
+c19 = newc;
+Circle(c19) = {p14, p10, p16};
+//+
+c20 = newc;
+Circle(c20) = {p16, p10, p13};
+//+
+
+
+//+
+cl3 = newcl;
+Curve Loop(cl3) = {c9,c17,c14};
+//+
+s2 = news;
+Surface(s2) = {cl3};
+//+
+cl4 = newcl;
+Curve Loop(cl4) = {c9,-c20,-c15};
+//+
+s3 = news;
+Surface(s3) = {cl4};
+//+
+cl5 = newcl;
+Curve Loop(cl5) = {c12,-c14,c18};
+//+
+s4 = news;
+Surface(s4) = {cl5};
+//+
+cl6 = newcl;
+Curve Loop(cl6) = {c12,c15,-c19};
+//+
+s5 = news;
+Surface(s5) = {cl6};
+//+
+cl7 = newcl;
+Curve Loop(cl7) = {c17,-c13,-c10};
+//+
+s6 = news;
+Surface(s6) = {cl7};
+//+
+cl8 = newcl;
+Curve Loop(cl8) = {c10,-c16,c20};
+//+
+s7 = news;
+Surface(s7) = {cl8};
+//+
+cl9 = newcl;
+Curve Loop(cl9) = {c13,c18,-c11};
+//+
+s8 = news;
+Surface(s8) = {cl9};
+//+
+cl10 = newcl;
+Curve Loop(cl10) = {-c11,-c16,-c19};
+//+
+s9 = news;
+Surface(s9) = {cl10};
+
+//+
+sl3 = newsl;
+Surface Loop(sl3) = {s2,s3,s4,s5,s6,s7,s8,s9};
+
+v4 = newv;
+Volume(v4) = {sl3,sl2,sl1};
+
+Physical Surface(3333) = {s2,s3,s4,s5,s6,s7,s8,s9};
+Physical Volume(3000) = v4;
+
+DefineConstant[
+  Flag_Infinity = {0, Choices{0,1},
+					Name "Input/01Use shell transformation to infinity"}
+];
+// radius for surrounding air with transformation to infinity
+If(Flag_Infinity==1)
+  label_Rext = "1Outer [m]";
+Else
+  label_Rext = "1[m]";
+EndIf
+
+Rint = 200*cm;
+Rext = 280*cm;
+
+DefineConstant[
+  Rint = {200*cm, Min 0.15, Max 0.9, Step 0.1, Name StrCat[pp2, "7Inner [m]"],
+    Visible (Flag_Infinity==1), Highlight Str[colorpp] },
+  Rext = {280*cm, Min Rint, Max 1, Step 0.1, Name StrCat[pp2, StrCat["8", label_Rext]],
+    Visible 1, Highlight Str[colorpp] }
+];
+/*
+lc0 = Pi*Rint/6; // air
+Characteristic Length { PointsOf{ Volume{v4}; } } = lc0;
+
+lc1  = 2*wcoil/2; // coil
+Characteristic Length { PointsOf{ Volume{e(1)}; } } = lc1;
+*/
+//Mesh 2;
+//Mesh 3;
+//Mesh.SaveAll = 1;
+//Save "mergedMesh.msh";

dG3D/benchmarks/weakCouplingSMP/run.py

@@ -32,8 +32,8 @@ finalDisp = 0.01;
 EMncycles = 1;
 
 # num of steps used in TM and EM solver
-TMnstep = 8;
-EMnstep = 80; #8*EMncycles;
+TMnstep = 16;
+EMnstep = 40; #8*EMncycles;
 
 stepIter = 0;
 
@@ -146,7 +146,7 @@ for currenttime in timeSteps:
   #fop.write('\nMerge "deformed_mesh_{}.msh"; '.format(TMnstep));
   fop.close();
 
-  os.system('gmsh merge_mesh.geo -3 -order 1 -o EMmergedMesh.msh')
+  os.system('gmsh new_merge_mesh.geo -3 -order 1 -o EMmergedMesh.msh')
   
   #if stepIter>1:
     # copy deformed mesh files to TM dir for backup

[Vinayak Gholap]

Note: Changed eqRatio in order to have a better convergence with abs residual converged value less than 1e-6.

[Vinayak Gholap]

Added new merge mesh file which:

• Takes new point, curve, curve loop, surface, surface loop and volume IDs automatically.
• Uses Transfinite Curve and Transfinite Surface for coil domain
• Extrudes coil domain using Layers

Earlier merge mesh file resulted in errors due to existing IDs for points and curves that were hand-coded numbers.