Fixes to code to work with elliptic problems again

examples/FEM4EllipticPDE_bhmag.cpp

@@ -149,6 +149,7 @@ int main(int argc, char**argv) {
         //Solver->SetSigmaFunction(implSigma);
         Solver->SetBoundaryStep(.05);
         Solver->SetGrid(MeshReader->GetOutput());
+        Solver->SetupPotential();
         //Solver->SetGrid(rGrid);
         Solver->Solve(argv[3]);
 

examples/borehole/sphere.geo

@@ -22,7 +22,7 @@ D0 = 10;          // Top of magnet
 D1 = 11;          // Bottom of magnet
 radius = 2.25;     // Radius of the damn thing
 
-lc = radius*10;   //  0.25;   // Target element size
+lc = radius;   //  0.25;   // Target element size
 
 // Total Solution Space
 Box = 10*radius; // The down side of potential
@@ -37,6 +37,31 @@ 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);
@@ -92,28 +117,6 @@ 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};}};
 
-/////////////////////////////////////
-// 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}; };
-
 /* This is GOOD */
 Surface{60} In Volume{v[1]};
 Surface{t1[0]} In Volume{v[1]};
@@ -128,13 +131,13 @@ 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[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("(.25 - F1)^2 + %g", cellSize );  // WORKS
-Field[2].F = Sprintf("(%g - F1)^2 + %g", radius, cellSize/2. );
-Background Field = 2;
+//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;

examples/borehole/sphere.sh

@@ -1,4 +1,9 @@
 
+gmsh  -3 -format vtk -o sphere.vtk  sphere.geo 
+gmsh  -2 -format stl -o sphereBox.stl  sphereBox.geo 
+
+
+
 #gmsh -3 -format vtk -o sphere.vtk  sphere.geo
 #../utVTKEdgeGsphere sphere.vtk .25 magnet  sphereG.vtu
 #../utFEM4EllipticPDE_bhmag  sphereG.vtu  sphereG.vtu   sphereOut.vtu 

examples/borehole/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 = 2.25;     // Radius of the damn thing
+
+lc = radius;   //  0.25;   // Target element size
+
+// Total Solution Space
+Box = 10*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

include/FEM4EllipticPDE.h

@@ -125,6 +125,10 @@ namespace Lemma {
              */
             void SetupDC(DCSurvey* Survey, const int& ij);
 
+            /** Sets up the potential problem from the VTK file
+             */
+            void SetupPotential();
+
             /** Performs solution */
             void Solve(const std::string& fname);
 

src/FEM4EllipticPDE.cpp

@@ -214,8 +214,9 @@ namespace Lemma {
 
         //Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
         cg.setMaxIterations(3000);
-        //cg.compute(A);
-        //std::cout << "Computed     " << std::endl;
+
+        std::cout << "A: "  << A.rows() << "\t" << A.cols() << std::endl;
+        std::cout << "g: "  << g.rows() << "\t" << g.cols() << std::endl;
 
         VectorXr u = cg.solve(g);
         std::cout << "#iterations:     " << cg.iterations() << std::endl;
@@ -257,20 +258,31 @@ namespace Lemma {
         std::cout << "Building Stiffness Matrix (A) " << std::endl;
         std::cout << "\tMesh has " << vtkGrid->GetNumberOfCells()  << " cells "  << std::endl;
 
-
   	    //Eigen::SparseMatrix<Real>
         A.resize(vtkGrid->GetNumberOfPoints(), vtkGrid->GetNumberOfPoints());
         std::vector< Eigen::Triplet<Real> > coeffs;
 
+        if ( !vtkGrid->GetPointData()->GetScalars("vtkValidPointMask") ) {
+            throw std::runtime_error("No vtkValidPointMask");
+        }
+
+        if ( !vtkGrid->GetCellData()->GetScalars("G") && !vtkGrid->GetPointData()->GetScalars("G") ) {
+            throw std::runtime_error("No Cell or Point Data G");
+        }
+
+        bool GCell = false;
+        if ( vtkGrid->GetCellData()->GetScalars("G") ) {
+            GCell = true;
+        }
+
+
         // Here we iterate over all of the cells
         for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
 
             assert ( vtkGrid->GetCell(ic)->GetNumberOfPoints() == 4 );
             // TODO, in production code we might not want to do this check here
             if ( vtkGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
-                std::cout << "DOOM FEM4EllipticPDE encountered non-tetrahedral mesh\n";
-                std::cout << "Number of points in cell " << vtkGrid->GetCell(ic)->GetNumberOfPoints() << std::endl ;
-                exit(1);
+                throw std::runtime_error("Non-tetrahedral mesh encountered!");
             }
 
             // construct coordinate matrix C
@@ -292,8 +304,16 @@ namespace Lemma {
                 ID[1] = Ids->GetId(1);
                 ID[2] = Ids->GetId(2);
                 ID[3] = Ids->GetId(3);
-            Real sum(0);
-            Real sigma_bar = vtkGrid->GetCellData()->GetScalars()->GetTuple1(ic);
+            Real sum(0), sigma_bar(0);
+            if (GCell) {
+                sigma_bar = vtkGrid->GetCellData()->GetScalars("G")->GetTuple1(ic);
+            } else {
+                sigma_bar  = vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0]);
+                sigma_bar += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1]);
+                sigma_bar += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2]);
+                sigma_bar += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3]);
+                sigma_bar /= 4.;
+            }
 
             for (int ip=0; ip<4; ++ip) {
                 for (int ip2=0; ip2<4; ++ip2) {
@@ -302,7 +322,6 @@ namespace Lemma {
                         //   solve for the boundaries? Is one better? This seems to work, which is nice.
                         //Real bdry = V*(1./(BndryH*BndryH))*BndrySigma*bndryCnt( ID[ip] ); // + sum;
                         Real bb = vtkGrid->GetPointData()->GetScalars("vtkValidPointMask")->GetTuple(ID[ip])[0];
-                        //std::cout << "bb " << bb  << std::endl;
                         Real bdry = V*(1./(BndryH*BndryH))*BndrySigma*bb; // + sum;
                         coeffs.push_back( Eigen::Triplet<Real> ( ID[ip], ID[ip2], bdry +  Phi.col(ip).tail<3>().dot(Phi.col(ip2).tail<3>() ) * V * sigma_bar ) );
                     } else {
@@ -315,6 +334,43 @@ namespace Lemma {
             std::cout <<  "\r" << (int)(1e2*((float)(ic) / (float)(vtkGrid->GetNumberOfCells()))) << std::flush ;
         }
         A.setFromTriplets(coeffs.begin(), coeffs.end());
+    }
+
+    void FEM4EllipticPDE::SetupPotential() {
+
+        ////////////////////////////////////////////////////////////
+	    // Load vector g
+        std::cout << "\nBuilding load vector (g)" << std::endl;
+        g = VectorXr::Zero(vtkGrid->GetNumberOfPoints());
+        std::cout << "made g with "  << vtkGrid->GetNumberOfPoints() << " cells" << std::endl;
+
+        for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
+
+                Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
+                for (int ip=0; ip<4; ++ip) {
+                    double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ip); //[ipc] ;
+                    C(ip, 0) = 1;
+                    C(ip, 1) =  pts[0];
+                    C(ip, 2) =  pts[1];
+                    C(ip, 3) =  pts[2];
+                }
+
+                Real V = (1./6.) * C.determinant();          // volume of tetrahedra
+                //Real V = C.determinant();          // volume of tetrahedra
+
+                vtkIdList* Ids = vtkGrid->GetCell(ic)->GetPointIds();
+                int ID[4];
+                    ID[0] = Ids->GetId(0);
+                    ID[1] = Ids->GetId(1);
+                    ID[2] = Ids->GetId(2);
+                    ID[3] = Ids->GetId(3);
+
+            for (int ip=0; ip<4; ++ip) {
+                 g(ID[ip]) +=  (V/4.) *  ( vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ip])[0] )  ;
+                 //if ( std::abs(vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ip])[0]) > 1e-3 )
+            }
+
+        }
 
     }