More work on the Linear Mag problem. Seems to be close.

examples/LinearMag/Sphere.cpp

@@ -31,6 +31,9 @@ int main( int argc, char** argv ) {
     std::cout << "# This program solves for the secondary magnetic  #\n";
     std::cout << "# field from a low susceptibily model in the      #\n";
     std::cout << "# presence of a homogeneous static inducing field #\n";
+    std::cout << "#                                                 #\n";
+    std::cout << "# (C) University of Utah, 2016                    #\n";
+    std::cout << "#                                                 #\n";
     std::cout << "###################################################\n\n";
 
     if (argc < 3) {
@@ -40,7 +43,7 @@ int main( int argc, char** argv ) {
     }
 
     LinearMag* Mag = LinearMag::New();
-        Mag->SetInducingMagField( 55234., 60, 14, NANOTESLA );
+        Mag->SetInducingMagField( 55234., 0, 0, NANOTESLA );
         Mag->SetVTUGridFile( argv[1] );
         Mag->CalculateRHS( "kappa" );
         Mag->Solve( argv[2] );

examples/LinearMag/sphere.geo

@@ -0,0 +1,148 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+radius = 3.25;     // Radius of the damn thing
+lc = radius/5;     //  0.25;   // Target element size
+
+// Total Solution Space
+Box = 3*radius; // The down side of potential
+X0 = -Box;
+X1 =  Box;
+Y0 = -Box;
+Y1 =  Box;
+Z0 = -Box;
+Z1 =  Box;
+
+cellSize=radius/10; ///10;
+dd = 0 ; //  1e-5; //cellSize; // .01;
+pio2=Pi/2;
+
+
+/////////////////////////////////////
+// Large Bounding box
+pp = newp;
+Point(pp)    = {X0, Y0, Z0, lc};
+Point(pp+1)  = {X1, Y0, Z0, lc};
+Point(pp+2)  = {X1, Y1, Z0, lc};
+Point(pp+3)  = {X0, Y1, Z0, lc};
+
+
+lv = newl;
+Line(lv) = {pp,pp+1};
+Line(lv+1) = {pp+1,pp+2};
+Line(lv+2) = {pp+2,pp+3};
+Line(lv+3) = {pp+3,pp};
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
+
+// Hard coded doom
+Plane Surface(125) = {lv+4};
+
+//v = newv;
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
+
+// Calculate offset effect
+theta = Asin(dd/radius);
+rr = radius * Cos(theta);
+
+///////////////////////////////////
+// Positive half sphere
+// create inner 1/8 shell
+p0 = newp;
+Point(p0)   = {      0,   0,       0, cellSize}; // origin
+Point(p0+1) = {    -rr,   0,      dd, cellSize};
+Point(p0+2) = {      0,  rr,      dd, cellSize};
+Point(p0+3) = {      0,   0,  radius, cellSize};
+Point(p0+4) = {      0,   0,      dd, cellSize}; // origin
+
+c0 = newc;
+Circle(c0  ) = {p0+1, p0+4, p0+2};       // Tricky, This one needs to be offset!
+Circle(c0+1) = {p0+3, p0, p0+1};
+Circle(c0+2) = {p0+3, p0, p0+2};
+
+Line Loop(10) = {c0, -(c0+2), c0+1} ;
+Ruled Surface (60) = {10};
+
+////////////////////////////////////////////////////////////
+// Negative half sphere
+p = newp;
+Point(  p) = {      0,      0,            0, cellSize};
+Point(p+1) = {    -rr,      0,          -dd, cellSize};
+Point(p+2) = {      0,     rr,          -dd, cellSize};
+Point(p+3) = {      0,      0,      -radius, cellSize};
+Point(p+4) = {      0,      0,          -dd, cellSize};
+
+cc = newc;
+Circle(cc  ) = {p+1, p+4, p+2};
+Circle(cc+1) = {p+3, p,   p+1};
+Circle(cc+2) = {p+3, p,   p+2};
+
+Circle(cc+3) = {p+3, p, p+2};
+Circle(cc+4) = {p+3, p, p+2};
+Circle(cc+5) = {p+3, p, p+2};
+
+ccl = newl;
+Line(ccl) = { p0+3, p+3 };
+
+Line Loop(11) = {cc, -(cc+2), cc+1} ;
+Ruled Surface (61) = {11};
+
+// create remaining 7/8 inner shells
+t1[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{60};}};
+t2[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{60};}};
+t3[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{60};}};
+//
+t4[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{61};}};
+t5[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{61};}};
+t6[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{61};}};
+
+/* This is GOOD */
+Surface{60} In Volume{v[1]};
+Surface{t1[0]} In Volume{v[1]};
+Surface{t2[0]} In Volume{v[1]};
+Surface{t3[0]} In Volume{v[1]};
+
+Surface{61} In Volume{v[1]};
+Surface{t4[0]} In Volume{v[1]};
+Surface{t5[0]} In Volume{v[1]};
+Surface{t6[0]} In Volume{v[1]};
+
+///////////////////////////////////////////////
+// Attractor Field
+
+// Field[1] = Attractor;
+// Field[1].NodesList = {p}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
+//
+// Field[2] = MathEval;
+// Field[2].F = Sprintf("(2.25 - F1)^1.01 + %g", radius/10 );  // WORKS
+//
+// Field[3] = MathEval;
+// Field[3].F = Sprintf("(12.25 - F1)^1.01 + %g", radius/5 );  // WORKS
+//
+// Field[4] = MathEval;
+// Field[4].F = Sprintf("(22.25 - F1)^1.01 + %g", radius );  // WORKS
+//
+// Field[5] = MathEval;
+// Field[5].F = Sprintf("(42.25 - F1)^1.01 + %g", radius*5 );  // WORKS
+//
+// //Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, 2*cellSize );
+// //Background Field = 2;
+//
+// // Finally, let's use the minimum of all the fields as the background mesh field
+// Field[7] = Min;
+// Field[7].FieldsList = {2, 3, 4, 5};
+// Background Field = 7;
+
+// Don't extend the elements sizes from the boundary inside the domain
+//Mesh.CharacteristicLengthExtendFromBoundary = 0;
+
+Physical Volume(1) = {v[1]};
+
+// To create the mesh run
+// gmsh sphere.gmsh -2 -v 0 -format msh -o sphere.msh
+//gmsh -3 -format msh1 -o outfile.msh sphere.geo

examples/LinearMag/sphere.sh

@@ -0,0 +1,8 @@
+#!/usr/bin/env bash
+gmsh  -3 -format vtk -o sphere.vtk  sphere.geo 
+gmsh  -2 -format stl -o sphereBox.stl  sphereBox.geo 
+../ResampleWithDataset  sphere.vtk  sphereBox.stl  MergedSphere.vtu 
+../VTKEdgeGsphere MergedSphere.vtu   3.25  mag  sphereMag.vtu
+./Sphere  sphereMag.vtu     sphereOutMag.vtu 
+
+#paraview --data=sphereOutGrav.vtu 

examples/LinearMag/sphereBox.geo

@@ -0,0 +1,136 @@
+/* This file is part of Lemma, a geophysical modelling and inversion API.
+ * More information is available at http://lemmasoftware.org
+ */
+
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+D0 = 10;          // Top of magnet
+D1 = 11;          // Bottom of magnet
+radius = 3.25;     // Radius of the damn thing
+
+lc = radius;   //  0.25;   // Target element size
+
+// Total Solution Space
+Box = 3*radius; // The down side of potential
+X0 = -Box;
+X1 =  Box;
+Y0 = -Box;
+Y1 =  Box;
+Z0 = -Box;
+Z1 =  Box;
+
+cellSize=lc; //300;
+dd = 0 ; //  1e-5; //cellSize; // .01;
+pio2=Pi/2;
+
+
+/////////////////////////////////////
+// Large Bounding box
+pp = newp;
+Point(pp)    = {X0, Y0, Z0, lc};
+Point(pp+1)  = {X1, Y0, Z0, lc};
+Point(pp+2)  = {X1, Y1, Z0, lc};
+Point(pp+3)  = {X0, Y1, Z0, lc};
+
+
+lv = newl;
+Line(lv) = {pp,pp+1};
+Line(lv+1) = {pp+1,pp+2};
+Line(lv+2) = {pp+2,pp+3};
+Line(lv+3) = {pp+3,pp};
+Line Loop(lv+4) = {lv, lv+1, lv+2, lv+3};
+
+// Hard coded doom
+Plane Surface(125) = {lv+4};
+
+//v = newv;
+v[] = Extrude {0, 0, Z1-Z0} { Surface{125}; };
+
+// // Calculate offset effect
+// theta = Asin(dd/radius);
+// rr = radius * Cos(theta);
+//
+// ///////////////////////////////////
+// // Positive half sphere
+// // create inner 1/8 shell
+// p0 = newp;
+// Point(p0)   = {      0,   0,       0, cellSize}; // origin
+// Point(p0+1) = {    -rr,   0,      dd, cellSize};
+// Point(p0+2) = {      0,  rr,      dd, cellSize};
+// Point(p0+3) = {      0,   0,  radius, cellSize};
+// Point(p0+4) = {      0,   0,      dd, cellSize}; // origin
+//
+// c0 = newc;
+// Circle(c0  ) = {p0+1, p0+4, p0+2};       // Tricky, This one needs to be offset!
+// Circle(c0+1) = {p0+3, p0, p0+1};
+// Circle(c0+2) = {p0+3, p0, p0+2};
+//
+// Line Loop(10) = {c0, -(c0+2), c0+1} ;
+// Ruled Surface (60) = {10};
+//
+// ////////////////////////////////////////////////////////////
+// // Negative half sphere
+// p = newp;
+// Point(  p) = {      0,      0,            0, cellSize};
+// Point(p+1) = {    -rr,      0,          -dd, cellSize};
+// Point(p+2) = {      0,     rr,          -dd, cellSize};
+// Point(p+3) = {      0,      0,      -radius, cellSize};
+// Point(p+4) = {      0,      0,          -dd, cellSize};
+//
+// cc = newc;
+// Circle(cc  ) = {p+1, p+4, p+2};
+// Circle(cc+1) = {p+3, p,   p+1};
+// Circle(cc+2) = {p+3, p,   p+2};
+//
+// Circle(cc+3) = {p+3, p, p+2};
+// Circle(cc+4) = {p+3, p, p+2};
+// Circle(cc+5) = {p+3, p, p+2};
+//
+// ccl = newl;
+// Line(ccl) = { p0+3, p+3 };
+//
+// Line Loop(11) = {cc, -(cc+2), cc+1} ;
+// Ruled Surface (61) = {11};
+//
+// // create remaining 7/8 inner shells
+// t1[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{60};}};
+// t2[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{60};}};
+// t3[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{60};}};
+// //
+// t4[] = Rotate {{0,0,1},{0,0,0},pio2  } {Duplicata{Surface{61};}};
+// t5[] = Rotate {{0,0,1},{0,0,0},pio2*2} {Duplicata{Surface{61};}};
+// t6[] = Rotate {{0,0,1},{0,0,0},pio2*3} {Duplicata{Surface{61};}};
+//
+// /* This is GOOD */
+// Surface{60} In Volume{v[1]};
+// Surface{t1[0]} In Volume{v[1]};
+// Surface{t2[0]} In Volume{v[1]};
+// Surface{t3[0]} In Volume{v[1]};
+//
+// Surface{61} In Volume{v[1]};
+// Surface{t4[0]} In Volume{v[1]};
+// Surface{t5[0]} In Volume{v[1]};
+// Surface{t6[0]} In Volume{v[1]};
+//
+// ///////////////////////////////////////////////
+// // Attractor Field
+//
+// //Field[1] = Attractor;
+// //Field[1].NodesList = {p}; //0, p0+1, p0+2, p0+3, p0+4, p, p+1, p+2, p+3, p+4};
+//
+// //Field[2] = MathEval;
+// //Field[2].F = Sprintf("(2.25 - F1)^2 + %g", cellSize*10 );  // WORKS
+// //Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, 2*cellSize );
+// //Background Field = 2;
+//
+// // Don't extend the elements sizes from the boundary inside the domain
+// //Mesh.CharacteristicLengthExtendFromBoundary = 0;
+//
+// Physical Volume(1) = {v[1]};
+
+// To create the mesh run
+// gmsh sphere.gmsh -2 -v 0 -format msh -o sphere.msh
+//gmsh -3 -format msh1 -o outfile.msh sphere.geo

examples/VTKEdgeGsphere.cpp

@@ -28,6 +28,8 @@
 #include <vtkXMLUnstructuredGridWriter.h>
 #include <vtkDoubleArray.h>
 #include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkCell.h>
 
 #include <cmath>
 #include <Eigen/Core>
@@ -42,7 +44,7 @@ int main(int argc, char**argv) {
         exit(EXIT_SUCCESS);
     }
 
-    vtkUnstructuredGrid* uGrid = vtkUnstructuredGrid::New();
+    vtkUnstructuredGrid* uGrid;// = vtkUnstructuredGrid::New();
 
     std::string fn = argv[1];
     if(fn.substr(fn.find_last_of(".") + 1) == "vtk") {
@@ -67,8 +69,8 @@ int main(int argc, char**argv) {
     }
 
 
-        int nc = uGrid->GetNumberOfCells()  ;
-        int nn = uGrid->GetNumberOfPoints()  ;
+    int nc = uGrid->GetNumberOfCells()  ;
+    int nn = uGrid->GetNumberOfPoints()  ;
 
     std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
 
@@ -84,80 +86,103 @@ int main(int argc, char**argv) {
         //G->SetName("G");
         if ( std::string(argv[3]) == "mag") {
             G->SetName("kappa");
+            G->SetNumberOfTuples(nc);
         }
         else {
             G->SetName("G");
+            G->SetNumberOfTuples(nn);
         }
     vtkDoubleArray* phi = vtkDoubleArray::New();
         phi->SetNumberOfComponents(1);
         phi->SetNumberOfTuples(nn);
         phi->SetName("analytic_phi");
 
-    // Loop over nodes, look at position, set with G if on boundary
     double point[3];
-    Eigen::Vector3d M; M << 0,0,1;
-    for (int in=0; in<nn; ++in) {
-        uGrid->GetPoint(in, &point[0]);
-        //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
-        double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
-
-        //double rho = sqrt( point[1]*point[1] + point[2]*point[2] );
-        // Calculate RHS
-        if ( std::string(argv[3]) == "gravity") {
-            if ( raddist < R + eps) {
-                G->InsertTuple1( in, 1 );
-            } else {
-                G->InsertTuple1( in, 0 );
+
+
+    if ( std::string(argv[3]) == "mag") {
+        // Loop over cells and decide if all of them lie within sphere
+        for (int ic=0; ic<nc; ++ic) {
+            assert ( uGrid->GetCell(ic)->GetNumberOfPoints() == 4 );
+            // TODO, in production code we might not want to do this check here
+            if ( uGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
+                throw std::runtime_error("Non-tetrahedral mesh encountered!");
             }
-        }
-        if ( std::string(argv[3]) == "magnet") {
-            if ( raddist > R - eps && raddist < R + eps ) {
-                // \rho = \nabla \cdot \mathbf{M}
-                // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
-                Eigen::Vector3d n;
-                n << point[0],point[1],point[2];
-                n.array() /= raddist;
-                G->InsertTuple1(in, n.dot(M) );
-            } else {
-                G->InsertTuple1( in, 0 );
+            // construct coordinate matrix C
+            //Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
+            bool cellin(true);
+            for (int ii=0; ii<4; ++ii) {
+                double* point =  uGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
+                if ( sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] ) > R+eps ) {
+                    cellin = false;
+                }
             }
-        }
-        if ( std::string(argv[3]) == "mag") {
-            if ( raddist < R + eps) {
-                G->InsertTuple1( in, 1 );
+            if (cellin) {
+                G->InsertTuple1( ic, 1 );
             } else {
-                G->InsertTuple1( in, 0 );
+                G->InsertTuple1( ic, 0 );
             }
         }
+        uGrid->GetCellData()->AddArray( G );
 
-        // Insert analytic phi part
-        /* magnetics problem p. 198 Jackson */
-        if (std::string(argv[3]) == "magnet") {
-            if (raddist < R) {
-                phi->InsertTuple1( in, (1./3.)*point[2] );
-            } else {
-                phi->InsertTuple1( in, 0);
-                double costheta = point[2]/raddist ;
-                //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) )  );
-                phi->InsertTuple1( in, (1./3.) * (R*R*R) * (costheta / (raddist*raddist))  );
+    } else {
+        // Loop over nodes, look at position, set with G if on boundary
+        Eigen::Vector3d M; M << 0,0,1;
+        for (int in=0; in<nn; ++in) {
+            uGrid->GetPoint(in, &point[0]);
+            //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
+            double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
+
+            //double rho = sqrt( point[1]*point[1] + point[2]*point[2] );
+            // Calculate RHS
+            if ( std::string(argv[3]) == "gravity") {
+                if ( raddist < R + eps) {
+                    G->InsertTuple1( in, 1 );
+                } else {
+                    G->InsertTuple1( in, 0 );
+                }
             }
-        } else if (std::string(argv[3]) == "gravity") {
-            double mass = 4./3. * PI * (R*R*R); // volume * density (1)
-            if (raddist < R) {
-                phi->InsertTuple1( in, mass * (( 3*(R*R) - raddist*raddist )/(2*(R*R*R))) ); // (1./3.)*point[2] );
-                //phi->InsertTuple1( in,  (2*PI/3)*(3*R*R-raddist*raddist)  ); // (1./3.)*point[2] );
-            } else {
-                //phi->InsertTuple1( in, 0);
-                //double costheta = point[2]/raddist ;
-                //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) )  );
-                phi->InsertTuple1( in, mass/raddist );
+            if ( std::string(argv[3]) == "magnet") {
+                if ( raddist > R - eps && raddist < R + eps ) {
+                    // \rho = \nabla \cdot \mathbf{M}
+                    // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
+                    Eigen::Vector3d n;
+                    n << point[0],point[1],point[2];
+                    n.array() /= raddist;
+                    G->InsertTuple1(in, n.dot(M) );
+                } else {
+                    G->InsertTuple1( in, 0 );
+                }
+            }
+            // Insert analytic phi part
+            /* magnetics problem p. 198 Jackson */
+            if (std::string(argv[3]) == "magnet") {
+                if (raddist < R) {
+                    phi->InsertTuple1( in, (1./3.)*point[2] );
+                } else {
+                    phi->InsertTuple1( in, 0);
+                    double costheta = point[2]/raddist ;
+                    //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) )  );
+                    phi->InsertTuple1( in, (1./3.) * (R*R*R) * (costheta / (raddist*raddist))  );
+                }
+            } else if (std::string(argv[3]) == "gravity") {
+                double mass = 4./3. * PI * (R*R*R); // volume * density (1)
+                if (raddist < R) {
+                    phi->InsertTuple1( in, mass * (( 3*(R*R) - raddist*raddist )/(2*(R*R*R))) ); // (1./3.)*point[2] );
+                    //phi->InsertTuple1( in,  (2*PI/3)*(3*R*R-raddist*raddist)  ); // (1./3.)*point[2] );
+                } else {
+                    //phi->InsertTuple1( in, 0);
+                    //double costheta = point[2]/raddist ;
+                    //phi->InsertTuple1( in, -(1./3.)*(R*R*R) * ( costheta / (rho*rho) )  );
+                    phi->InsertTuple1( in, mass/raddist );
+                }
             }
         }
-    }
 
-    // Add new data
-    uGrid->GetPointData()->AddArray( phi );
-    uGrid->GetPointData()->SetScalars( G );
+        // Add new data
+        uGrid->GetPointData()->AddArray( phi );
+        uGrid->GetPointData()->SetScalars( G );
+    }
 
     // Write out new file
     vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();

src/FEM4EllipticPDE.cpp

@@ -145,11 +145,13 @@ namespace Lemma {
     //      Method:  SetVTUGridFile
     //--------------------------------------------------------------------------------------
     void FEM4EllipticPDE::SetVTUGridFile ( std::string& fname  ) {
+        std::cout  << "Loading VTK .vtu file " << fname << std::endl;
         vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
         MeshReader->SetFileName(fname.c_str());
         MeshReader->Update();
         //vtkGrid->DeepCopy( MeshReader->GetOutput() );
-        vtkGrid->ShallowCopy( MeshReader->GetOutput() );
+        //vtkGrid->ShallowCopy( MeshReader->GetOutput() );
+        vtkGrid = MeshReader->GetOutput();
         MeshReader->Delete();
         return ;
     }		// -----  end of method FEM4EllipticPDE::SetVTKGridFile  -----
@@ -164,8 +166,8 @@ namespace Lemma {
         MeshReader->SetFileName(fname);
         MeshReader->Update();
         //vtkGrid->DeepCopy( MeshReader->GetOutput() );
-        vtkGrid = MeshReader->GetOutput();
         //vtkGrid->ShallowCopy( MeshReader->GetOutput() );
+        vtkGrid = MeshReader->GetOutput();
         MeshReader->Delete();
         return ;
     }		// -----  end of method FEM4EllipticPDE::SetVTKGridFile  -----
@@ -547,10 +549,8 @@ namespace Lemma {
 	    // Load vector g
         std::cout << "\nBuilding load vector (g)" << std::endl;
         g = VectorXr::Zero(vtkGrid->GetNumberOfPoints());
-        std::cout << "made g with "  << vtkGrid->GetNumberOfPoints() << " points" << std::endl;
 
         for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
-
             Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
             for (int ii=0; ii<4; ++ii) {
                 double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
@@ -573,9 +573,9 @@ namespace Lemma {
 
             // Apply Dirichlet conditions
             for (int ii=0; ii<4; ++ii) {
-                if (vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0]) {
-                    g(ID[ii]) += vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
-                }
+                //if (vtkGrid->GetPointData()->GetScalars("InHomogeneousDirichlet")->GetTuple(ID[ii])[0]) {
+                //    g(ID[ii]) += vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
+                //}
             }
 
             /* the 4 faces of the tetrahedra
@@ -632,6 +632,7 @@ namespace Lemma {
                 g(ID[3]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])*TA/3. * i4pi ;
             }
         }
+        std::cout << "made g with "  << vtkGrid->GetNumberOfPoints() << " points" << std::endl;
     }
 
     #if 0

src/LinearMag.cpp

@@ -173,14 +173,27 @@ void LinearMag::ScaleB0 ( const MAGUNITS& U ) {
 //      Method:  CalculateRHS
 //--------------------------------------------------------------------------------------
 void LinearMag::CalculateRHS ( const std::string& susName ) {
+
+    if ( !vtkGrid->GetCellData()->GetScalars(susName.c_str()) ) {
+        std::string err("No cell data by name ");
+        err.append(susName);
+        throw std::runtime_error(err.c_str());
+    }
+
     vtkDoubleArray* G = vtkDoubleArray::New();
     G->SetNumberOfComponents(1);
     G->SetNumberOfTuples( vtkGrid->GetNumberOfPoints() );
     G->SetName("G");
+    //g.resize(vtkGrid->GetNumberOfPoints());
+    VectorXr GG = VectorXr::Zero( vtkGrid->GetNumberOfPoints() );
 
     // Iterate over all the points or all of the cells?
     for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
 
+        if ( vtkGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
+            throw std::runtime_error("Non-tetrahedral mesh encountered!");
+        }
+
         Real cellSus = vtkGrid->GetCellData()->GetScalars(susName.c_str())->GetTuple(ic)[0];
 
         Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
@@ -207,23 +220,59 @@ void LinearMag::CalculateRHS ( const std::string& susName ) {
             ID[1] ID[2] ID[3]
         */
         // Face 0, ID 0,1,2
-        /*
         Eigen::Matrix<Real, 3, 2> CC = Eigen::Matrix<Real, 3, 2>::Ones() ;
         {
+
             CC.col(1) = C.row(0).tail<3>() - C.row(1).tail<3>();
             CC.col(2) = C.row(0).tail<3>() - C.row(2).tail<3>();
             Vector3r nhat = CC.col(1).cross(CC.col(2));
             nhat.array() /= nhat.norm();
             Real flux = cellSus*nhat.dot(B0);
-            g(ID[0]) += flux;
-            g(ID[1]) += flux;
-            g(ID[2]) += flux;
-                    //  vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])*TA/3. * i4pi ;
+            GG(ID[0]) += flux;
+            GG(ID[1]) += flux;
+            GG(ID[2]) += flux;
         }
-        */
+        // Face 1, ID 0,1,3
+        {
+            CC.col(1) = C.row(0).tail<3>() - C.row(1).tail<3>();
+            CC.col(2) = C.row(0).tail<3>() - C.row(3).tail<3>();
+            Vector3r nhat = CC.col(1).cross(CC.col(2));
+            nhat.array() /= nhat.norm();
+            Real flux = cellSus*nhat.dot(B0);
+            GG(ID[0]) += flux;
+            GG(ID[1]) += flux;
+            GG(ID[3]) += flux;
+        }
+        // Face 1, ID 0,2,3
+        {
+            CC.col(1) = C.row(0).tail<3>() - C.row(2).tail<3>();
+            CC.col(2) = C.row(0).tail<3>() - C.row(3).tail<3>();
+            Vector3r nhat = CC.col(1).cross(CC.col(2));
+            nhat.array() /= nhat.norm();
+            Real flux = cellSus*nhat.dot(B0);
+            GG(ID[0]) += flux;
+            GG(ID[2]) += flux;
+            GG(ID[3]) += flux;
+        }
+        // Face 1, ID 1,2,3
+        {
+            CC.col(1) = C.row(1).tail<3>() - C.row(2).tail<3>();
+            CC.col(2) = C.row(1).tail<3>() - C.row(3).tail<3>();
+            Vector3r nhat = CC.col(1).cross(CC.col(2));
+            nhat.array() /= nhat.norm();
+            Real flux = cellSus*nhat.dot(B0);
+            GG(ID[1]) += flux;
+            GG(ID[2]) += flux;
+            GG(ID[3]) += flux;
+        }
+
+    }
+    for (int ip=0; ip<vtkGrid->GetNumberOfPoints(); ++ip) {
+        G->InsertTuple1( ip, GG(ip) );
     }
 
     vtkGrid->GetPointData()->AddArray( G );
+    vtkGrid->GetPointData()->SetScalars( G );
     return ;
 }		// -----  end of method LinearMag::CalculateRHS  -----