nothing wrong with non-div-free RHS... same story as non-div-free guess from a...

nothing wrong with non-div-free RHS... same story as non-div-free guess from a gmres for instance...

example/maxwell_sibc_divfree/Makefile

-BIN=../../bin
-
-sphere: init cim
-
-init:
-	gmsh sphere.geo -3
-
-cim:
-	python ${BIN}/cim.py sphere.pro sphere.msh Solve 8.22543e9+2.46361e1j 1e8 -nodes 10 -lStart 4
-
-clean:
-	rm -f *.pos
-	rm -f *.msh
-	rm -f *.pre
-	rm -f *.res

example/maxwell_sibc_divfree/README.md

-Same example as [maxwell\_sibc](../maxwell_sibc), but imposes a divergence free
-source for Beyn. This is achieved by taking the curl of the random Beyn source.
-
-This example requires a [GetDP](http://getdp.info) version compiled with [Gmsh](
-http://gmsh.info) to handle Fields.

example/maxwell_sibc_divfree/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[x() = {},  // Solution
-                 b() = {}]; // Right hand side
-
-  // Control data //
-  DefineConstant[imposeRHS = 0,          // Should I use an imposed RHS?
-                 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_divfree/cimResolution.pro

-/////////////////////////////////////////////////////////////////////
-// Resolution for contour integral method.                         //
-// Transforms Beyn random source into a divergence free one.       //
-// 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; }
-      { Name B; NameOfFormulation Source;         Type ComplexValue; }
-    }
-    Operation{
-      // Divergence free source
-      If(imposeRHS)
-        Generate[B];
-        CopyRightHandSide[b(), B];
-        Solve[B];
-        PostOperation[Source];
-      EndIf
-
-      // Beyn (source given by dSrc)
-      Generate[A];
-
-      If(!doPostpro && !doApply)
-        Solve[A];
-        CopySolution[A, x()];
-      EndIf
-
-      If(doApply)
-        CopySolution[x(), A];
-        Apply[A];
-        CopySolution[A, x()];
-      EndIf
-
-      If(doPostpro)
-        CopySolution[x(), A];
-        PostOperation[Post];
-      EndIf
-    }
-  }
-}

example/maxwell_sibc_divfree/sphere.dat

-// Geometry //
-DefineConstant[R  = { 100e-3, Name "Input/00Geometry/00Radius [m]" }];
-DefineConstant[K1 = { 1/R,    Name "Input/00Geometry/01Curvature 1 [m^-1]",
-                              ReadOnly 1}];
-DefineConstant[K2 = { 1/R,    Name "Input/00Geometry/02Curvature 2 [m^-1]",
-                              ReadOnly 1}];
-
-// Mesh //
-DefineConstant[md = { 10,     Name "Input/01Mesh/00Density [Radius^-1]" }];
-cl = R/md;
-
-// Material //
-DefineConstant[Sigma = { 1e15,
-                         Name "Input/02Material/00Conductivity [S m^-1]" }];

example/maxwell_sibc_divfree/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_divfree/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) //
-  Sigma[] = Sigma; // [S/m]
-
-  // Leontovich SIBC (H-Formulation) //
-  SL[] = Sqrt[(J[] * Re[Puls[]] * Mu0) / Sigma[]]; // H-Formulation
-  SIbc[] = -1.0 / SL[];                            // E-Formulation
-
-  // Curl of Beyn source //
-  dSrc[] = ComplexVectorField[XYZ[]]{1};
-}
-
-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 Source; Type FemEquation;
-    Quantity{
-      { Name s; Type Local; NameOfSpace HCurl; }
-    }
-    Equation{
-      Galerkin{ [Dof{s}, {s}];
-        In Omega; Jacobian Jac; Integration I1; }
-    }
-  }
-  { 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; }
-
-      // Source
-      If(imposeRHS)
-        Galerkin{ [dSrc[], {e}];
-          In Omega; Jacobian Jac; Integration I1; }
-      EndIf
-    }
-  }
-}
-
-Include "./cimResolution.pro"; // Resolution for cim.py
-
-PostProcessing{
-  { Name Source; NameOfFormulation Source;
-    Quantity{
-      { Name dSrc; Value{ Term{ [{d s}]; In Omega; Jacobian Jac; } } }
-    }
-  }
-  { Name Post;   NameOfFormulation FormulationCIM;
-    Quantity{
-      { Name E;    Value{ Term{ [{  e}]; In Omega; Jacobian Jac; } } }
-    }
-  }
-}
-
-PostOperation{
-  { Name Source; NameOfPostProcessing Source;
-    Operation{
-      Print[dSrc, OnElementsOf Omega, StoreInField 1];
-    }
-  }
-  { Name Post;   NameOfPostProcessing Post;
-    Operation{
-      Print[E,    OnElementsOf Omega, File fileName];
-    }
-  }
-}