new test case with lagrange correction

example/maxwell_sibc_lagrange/Makefile

+BIN=../../bin
+
+sphere: init cim
+
+init:
+	gmsh sphere.geo -3
+
+cim:
+	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1e10+1e9j 4e9 -nodes 20 -lStart 20 -setnumber isSC 0
+
+#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 7.4e9+7.3e8j 1e9 -nodes 10 -lStart 4 -setnumber isSC 0
+#7.474355560164e+09 + 7.742536715699e+08j
+
+#	python ${BIN}/cim.py sphere.pro sphere.msh Solve 1.0e10+1.1e9j 1e9 -nodes 10 -lStart 6 -setnumber isSC 0
+#1.050336615299e+10 + 1.186115651023e+09j
+
+clean:
+	rm -f *.pos
+	rm -f *.msh
+	rm -f *.pre
+	rm -f *.res

example/maxwell_sibc_lagrange/cimParameters.pro

+/////////////////////////////////////////////////////////////////
+// Parameters for contour integral method.                     //
+// The master file should inclue "cimResolution.pro" later on. //
+/////////////////////////////////////////////////////////////////
+Function{
+  // Physical data //
+  DefineConstant[angularFreqRe = 1,  // [rad/s]
+                 angularFreqIm = 0]; // [rad/s]
+
+  Puls[]  = Complex[angularFreqRe,
+                    angularFreqIm];  // [rad/m]
+
+  // Algebraic data //
+  DefineConstant[nRHS = 1];       // Number of RHS for this run
+  For i In {0:nRHS-1}
+    DefineConstant[x~{i}() = {},  // Solution
+                   b~{i}() = {}]; // Right hand side
+  EndFor
+
+  // Control data //
+  DefineConstant[doInit    = 0,          // Should I initialize some stuff?
+                 doSolve   = 0,          // Should I solve Ax = b?
+                 doPostpro = 0,          // Should I only create a view for x()?
+                 doApply   = 0,          // Should I only apply x(): x <- Ax?
+                 fileName  = "eig.pos"]; // Postpro file name
+}

example/maxwell_sibc_lagrange/cimResolution.pro

+/////////////////////////////////////////////////////////////////////
+// Resolution for contour integral method.                         //
+// Warning, the master file should include "cimParameters.pro",    //
+// and should also define "FormulationCIM" the formulation to use. //
+/////////////////////////////////////////////////////////////////////
+Resolution{
+  { Name Solve;
+    System{ { Name A; NameOfFormulation FormulationCIM; Type ComplexValue; } }
+    Operation{
+      Generate[A];
+
+      If(doInit)
+        CopySolution[A, x~{0}()];
+      EndIf
+
+      If(doSolve)
+        // Full solve for first RHS
+        CopyRightHandSide[b~{0}(), A];
+        Solve[A];
+        CopySolution[A, x~{0}()];
+
+        // SolveAgain for other RHSs
+        For i In {1:nRHS-1}
+          CopyRightHandSide[b~{i}(), A];
+          SolveAgain[A];
+          CopySolution[A, x~{i}()];
+        EndFor
+      EndIf
+
+      If(doApply)
+        CopySolution[x~{0}(), A];
+        Apply[A];
+        CopySolution[A, x~{0}()];
+      EndIf
+
+      If(doPostpro)
+        CopySolution[x~{0}(), A];
+        PostOperation[Post];
+      EndIf
+    }
+  }
+}

example/maxwell_sibc_lagrange/sphere.dat

+// Geometry //
+DefineConstant[R  = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
+
+// Mesh //
+DefineConstant[md = { 10,     Name "Input/01Mesh/00Density [Radius^-1]" }];
+cl = R/md;
+
+// Material //
+DefineConstant[isSC   = { 0,
+                          Name "Input/02Material/00Is superconductor?",
+                          GmshOption "Reset",
+                          Choices {0 = "No", 1 = "Yes"} }];
+If(isSC == 0)
+DefineConstant[Sigma  = { 1e00,
+                          Name "Input/02Material/01Conductivity [S m^-1]" }];
+EndIf
+If(isSC == 1)
+DefineConstant[Sigma  = { 1e9,
+                          Name "Input/02Material/01Conductivity [S m^-1]" }];
+DefineConstant[Lambda = { 1e-8,
+                          Name "Input/02Material/02London depth [m]"      }];
+EndIf

example/maxwell_sibc_lagrange/sphere.geo

+Include "sphere.dat";
+
+Point(0) = { 0,  0,  0, cl};
+
+Point(1) = {+R,  0,  0, cl};
+Point(2) = { 0, +R,  0, cl};
+Point(3) = {-R,  0,  0, cl};
+Point(4) = { 0, -R,  0, cl};
+Point(5) = { 0,  0, +R, cl};
+Point(6) = { 0,  0, -R, cl};
+
+Circle(01) = {1, 0, 2};
+Circle(02) = {2, 0, 3};
+Circle(03) = {3, 0, 4};
+Circle(04) = {4, 0, 1};
+
+Circle(05) = {5, 0, 1};
+Circle(06) = {1, 0, 6};
+Circle(07) = {6, 0, 3};
+Circle(08) = {3, 0, 5};
+
+Circle(09) = {2, 0, 6};
+Circle(10) = {6, 0, 4};
+Circle(11) = {4, 0, 5};
+Circle(12) = {5, 0, 2};
+
+Line Loop(1) = {1, -12, 5};
+Line Loop(2) = {1, 9, -6};
+Line Loop(3) = {9, 7, -2};
+Line Loop(4) = {2, 8, 12};
+Line Loop(5) = {11, 5, -4};
+Line Loop(6) = {4, 6, 10};
+Line Loop(7) = {10, -3, -7};
+Line Loop(8) = {3, 11, -8};
+
+Ruled Surface(1) = {1};
+Ruled Surface(2) = {2};
+Ruled Surface(3) = {3};
+Ruled Surface(4) = {4};
+Ruled Surface(5) = {5};
+Ruled Surface(6) = {6};
+Ruled Surface(7) = {7};
+Ruled Surface(8) = {8};
+
+Surface Loop(1) = {3, 2, 1, 4, 8, 7, 6, 5};
+Volume(1) = {1};
+
+Physical  Volume(1) = {1};
+Physical Surface(2) = {1, 2, 3, 4, 5, 6, 7, 8};
+Physical   Point(3) = {0};

example/maxwell_sibc_lagrange/sphere.pro

+Include "sphere.dat";
+
+Group{
+  Omega = Region[1];
+  Bndry = Region[2];
+}
+
+Include "./cimParameters.pro"; // Functions for cim.py
+
+Function{
+  // Imaginary Unit //
+  J[] = Complex[0, 1];
+
+  // Physical Constants //
+  C0   = 299792458;      // [m/s]
+  Mu0  = 4*Pi*1e-7;      // [H/m]
+  Eps0 = 1/(Mu0 * C0^2); // [F/m]
+
+  // Conductivity (coming from sphere.dat) //
+  If(isSC == 1)                                            // Superconductor
+    Sigma[] = Sigma - J[] / (Mu0 * Lambda^2 * Re[Puls[]]); // [S/m]
+  EndIf
+  If(isSC == 0)                                            // Normal conductor
+    Sigma[] = Sigma;                                       // [S/m]
+  EndIf
+
+  // Leontovich SIBC (H-Formulation) //
+  PulsG1[]   = 7.474355560164e+09 + 7.742536715699e+08*J[];
+  PulsG2[]   = 1.050336615299e+10 + 1.186115651023e+09*J[];
+
+  PulsBar1[] = Re[PulsG1[]] - J[]*Im[PulsG1[]];
+  PulsBar2[] = Re[PulsG2[]] - J[]*Im[PulsG2[]];
+
+  L1[] = (Puls[] - PulsG2[]) / (PulsG1[] - PulsG2[]);
+  L2[] = (Puls[] - PulsG1[]) / (PulsG2[] - PulsG1[]);
+
+  //PulsCor[] = PulsBar1[]*Puls[]/PulsG1[]; //Re[Puls[]] + J[]*Im[Puls[]];
+  PulsCor[] = PulsBar1[]*L1[] + PulsBar2[]*L2[];
+  SL[] = Sqrt[(J[] * (Puls[]+PulsCor[])/2 * Mu0) / Sigma[]]; // H-Formulation
+  SIbc[] = -1.0 / SL[];                            // E-Formulation
+}
+
+Integration{
+  { Name I1;
+    Case{
+      { Type Gauss;
+        Case{
+          { GeoElement Tetrahedron; NumberOfPoints 4; } // O1 * O1
+          { GeoElement Triangle;    NumberOfPoints 3; } // O1 * O1
+          { GeoElement Line;        NumberOfPoints 2; } // O1 * O1
+        }
+      }
+    }
+  }
+}
+
+Jacobian{
+  { Name Jac;
+    Case{
+      { Region Omega; Jacobian Vol; }
+      { Region Bndry; Jacobian Sur; }
+    }
+  }
+}
+
+FunctionSpace{
+  { Name HCurl; Type Form1;
+    BasisFunction{
+      { Name sx; NameOfCoef ux; Function BF_Edge;
+        Support Region[{Omega, Bndry}]; Entity EdgesOf[All]; }
+    }
+  }
+}
+
+Formulation{
+  { Name FormulationCIM; Type FemEquation;
+    Quantity{
+      { Name e; Type Local; NameOfSpace HCurl; }
+    }
+    Equation{
+      // Maxwell
+      Galerkin{ [          1/Mu0 * Dof{d e}, {d e}];
+        In Omega; Jacobian Jac; Integration I1; }
+      Galerkin{ [      -Puls[]^2 * Eps0 * Dof{  e}, {  e}];
+        In Omega; Jacobian Jac; Integration I1; }
+
+      // IBC
+      Galerkin{ [-J[] * Puls[]   * SIbc[] * Dof{  e}, {  e}];
+        In Bndry; Jacobian Jac; Integration I1; }
+    }
+  }
+}
+
+Include "./cimResolution.pro"; // Resolution for cim.py
+
+PostProcessing{
+  { Name Post; NameOfFormulation FormulationCIM;
+    Quantity{
+      { Name E; Value{ Term{ [{e}]; In Omega; Jacobian Jac; } } }
+    }
+  }
+}
+
+PostOperation{
+  { Name Post; NameOfPostProcessing Post;
+    Operation{
+      Print[E, OnElementsOf Omega, File fileName];
+    }
+  }
+}