Cleaning up code, gravity seems to be working.

examples/VTKGsphere.cpp

@@ -113,15 +113,16 @@ int main(int argc, char**argv) {
                 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)
+            double mass = (4.*PI/3.) * (R*R*R); // volume * density (1)
+            double G = 1/(4.*PI);
             if (raddist < R) {
-                phi->InsertTuple1( in, mass * (( 3*(R*R) - raddist*raddist )/(2*(R*R*R))) ); // (1./3.)*point[2] );
+                phi->InsertTuple1( in, mass * G * (( 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 );
+                phi->InsertTuple1( in, G * mass/raddist );
             }
         }
     }

examples/borehole/sphere.geo

@@ -7,11 +7,11 @@
  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  */
 
-radius = 4.25;     // Radius of the damn thing
-lc = radius/5;     //  0.25;   // Target element size
+radius = 2.25;     // Radius of the damn thing
+lc = radius/6;     //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 5*radius; // The down side of potential
+Box = 3*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;
@@ -19,7 +19,7 @@ Y1 =  Box;
 Z0 = -Box;
 Z1 =  Box;
 
-cellSize=radius/10; ///10;
+cellSize=radius/6; ///10;
 dd = 0 ; //  1e-5; //cellSize; // .01;
 pio2=Pi/2;
 

examples/borehole/sphere.sh

@@ -2,7 +2,7 @@
 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 
-../VTKGsphere MergedSphere.vtu   4.25  gravity  sphereGrav.vtu
+../VTKGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
 ../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
 
 #paraview --data=sphereOutGrav.vtu 

examples/borehole/sphereBox.geo

@@ -9,12 +9,12 @@
 
 D0 = 10;          // Top of magnet
 D1 = 11;          // Bottom of magnet
-radius = 4.25;     // Radius of the damn thing
+radius = 2.25;     // Radius of the damn thing
 
 lc = radius;   //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 5*radius; // The down side of potential
+Box = 3*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;

include/FEM4EllipticPDE.h

@@ -134,7 +134,7 @@ namespace Lemma {
             void Solve(const std::string& fname);
 
             /** Performs solution */
-            void SolveOLD(const std::string& fname);
+            void SolveOLD2(const std::string& fname);
 
             // ====================  ACCESS        =======================
 

src/FEM4EllipticPDE.cpp

@@ -205,30 +205,38 @@ namespace Lemma {
         //ConstructLoadVector();
 
         std::cout << "\nSolving" << std::endl;
-        std::cout << "   rows\tcols\n";
-        std::cout << "A: "  << A.rows() << "\t" << A.cols() << std::endl;
-        std::cout << "g: "  << g.rows() << "\t" << g.cols() << std::endl;
+
+        std::cout << std::setw(5) << "  " << std::setw(14) << "rows"     << std::setw(14) << "cols"     << std::endl;
+        std::cout << std::setw(5) << "  " << std::setw(14) << "--------" << std::setw(14) << "--------" << std::endl;
+        std::cout << std::setw(5) << "A:" << std::setw(14) << A.rows()   << std::setw(14) << A.cols()   << std::endl;
+        std::cout << std::setw(5) << "g:" << std::setw(14) << g.rows()   << std::setw(14) << g.cols()   << std::endl;
 
         ////////////////////////////////////////////////////////////
         // Solving:
         //Eigen::SimplicialCholesky<Eigen::SparseMatrix<Real>, Eigen::Lower > chol(A);  // performs a Cholesky factorization of A
         //VectorXr u = chol.solve(g);
-
-        //Eigen::SparseLU<Eigen::SparseMatrix<Real, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > solver;
-        //solver.analyzePattern(A);
-        //solver.factorize(A);
-        //VectorXr u = solver.solve(g);
-
-        //#ifdef CGSOLVE
-        //Eigen::ConjugateGradient<Eigen::SparseMatrix<Real, Eigen::Lower >  Eigen::DiagonalPreconditioner > cg;
-        //Eigen::ConjugateGradient< Eigen::SparseMatrix<Real>  > cg(A);
+        //#define LUSOLVE
+        #ifdef LUSOLVE
+        Eigen::SparseLU<Eigen::SparseMatrix<Real, Eigen::ColMajor>, Eigen::COLAMDOrdering<int> > solver;
+        std::cout << "LU analyze pattern" << std::endl;
+        solver.analyzePattern(A);
+        std::cout << "LU factorizing" << std::endl;
+        solver.factorize(A);
+        std::cout << "LU solving" << std::endl;
+        solver.factorize(A);
+        VectorXr u = solver.solve(g);
+        #endif
+
+        #define CGSOLVE
+        #ifdef CGSOLVE
+        // TODO try IncompleteLUT preconditioner
         Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
         //cg.setMaxIterations(3000);
         //cg.setTolerance(1e-28);
         VectorXr u = cg.solve(g);
         std::cout << "#iterations:     " << cg.iterations() << std::endl;
         std::cout << "estimated error: " << cg.error()      << std::endl;
-        //#endif
+        #endif
 
         vtkDoubleArray *gArray = vtkDoubleArray::New();
         vtkDoubleArray *uArray = vtkDoubleArray::New();
@@ -284,7 +292,6 @@ namespace Lemma {
             GCell = true;
         }
 
-
         // Here we iterate over all of the cells
         for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
 
@@ -315,19 +322,8 @@ namespace Lemma {
                 ID[3] = Ids->GetId(3);
 
             Real 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.;
-            }
-            */
-
             // TEST VOID IN SIGMA!! TODO DON"T KEEP THIS
+            /*
             Real xc = C.col(1).array().mean();
             Real yc = C.col(2).array().mean();
             Real zc = C.col(3).array().mean();
@@ -336,40 +332,26 @@ namespace Lemma {
             } else {
                 sigma_bar = 1.;
             }
+            */
             sigma_bar = 1.;
 
 
             for (int ii=0; ii<4; ++ii) {
                 int bbi = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0];
                 if (bbi) {
+                    /* Dirichlet boundary */
                     coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[ii], 1));
                 } else {
                     for (int jj=0; jj<4; ++jj) {
-                        /* homogeneous Dirichlet boundary */
-                        //if (jj==ii && C(ii,3)==13.5) {
-                        //if (jj==ii && ((bb-1.)<1e-8) ) {
-                            // Apply Homogeneous Dirichlet Boundary conditions
-                        //    Real bb = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0];
-                        //    Real bdry = (1./(BndryH*BndryH))*BndrySigma*bb; // + sum;
-                            //Real bdry = GradPhi.col(ii).tail<3>().dot(GradPhi.col(ii).tail<3>())*BndrySigma*bb; // + sum;
-                        //    coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], bdry + GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );
-                            //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], 1));
-                            //break;
-                        //} else {
                         coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj],        GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );
-                        //}
-                    /* */
-                    /* Inhomogeneous Dirichlet */
-                    //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>())*V*sigma_bar ));
-                    // Stiffness matrix no longer contains boundary conditions...
-                    //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], GradPhi.col(ii).tail<3>().dot(GradPhi.col(jj).tail<3>() ) * V * sigma_bar ) );
+                    }
                 }
             }
-            }
-            //std::cout <<  "\r" << (int)(1e2*((float)(ic) / (float)(vtkGrid->GetNumberOfCells()))) << std::flush ;
+            std::cout <<  "\r" << (int)(1e2*((float)(ic) / (float)(vtkGrid->GetNumberOfCells()))) << std::flush ;
         }
         A.setFromTriplets(coeffs.begin(), coeffs.end());
-        //std::cout << "A\n" << A << std::endl;
+        A.finalize();
+        A.makeCompressed();
     }
 
     void FEM4EllipticPDE::SetupPotential() {
@@ -379,7 +361,7 @@ namespace Lemma {
         std::cout << "\nBuilding load vector (g)" << std::endl;
         g = VectorXr::Zero(vtkGrid->GetNumberOfPoints());
         std::cout << "made g with "  << vtkGrid->GetNumberOfPoints() << " points" << std::endl;
-        VectorXr DB = VectorXr::Zero(vtkGrid->GetNumberOfPoints());
+
         for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
 
             Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
@@ -392,8 +374,9 @@ namespace Lemma {
             }
 
             Real V = (1./6.) * C.determinant();              // volume of tetrahedra
-            Eigen::Matrix<Real, 4, 4> GradPhi = C.inverse(); // nabla \phi
+            //Eigen::Matrix<Real, 4, 4> GradPhi = C.inverse(); // nabla \phi
 
+            /* The indices */
             vtkIdList* Ids = vtkGrid->GetCell(ic)->GetPointIds();
             int ID[4];
                 ID[0] = Ids->GetId(0);
@@ -401,56 +384,19 @@ namespace Lemma {
                 ID[2] = Ids->GetId(2);
                 ID[3] = Ids->GetId(3);
 
-            /*
-            Real avg(0);
-            Real GG[4];
-            for (int ii=0; ii<4; ++ii) {
-                GG[ii] = vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0];
-                //avg /= 4.;
-            }
-            if ( std::abs( (GG[0]+GG[1]+GG[2]+GG[3])/4. - GG[0])  < 1e-5) {
-                avg = GG[0];
-            }
-            */
-
-            /*
-            VectorXr W = VectorXr::Zero(4);
+            /* Fill the RHS vector with Dirichlet conditions or fuction value */
             for (int ii=0; ii<4; ++ii) {
-                W[ii] = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0] *
-                        vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
-                DB[ID[ii]] = vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0] *
-                             vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
-            }
-            //auto G = GradPhi.block<3,4>(1,0).transpose()*GradPhi.block<3,4>(1,0)*W;
-            VectorXr G = GradPhi.block<3,4>(1,0).transpose()*GradPhi.block<3,4>(1,0)*W;
-            Real  sigma_bar(1.);
-            */
-
-            for (int ii=0; ii<4; ++ii) {
-            //     avg += vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0];
-            //     //if ( std::abs(vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0]) > 1e-3 )
-            //}
-            //TODO this seems wrong!
-            //avg /= 4.;
-                // TODO test code remove
                 if (vtkGrid->GetPointData()->GetScalars("HomogeneousDirichlet")->GetTuple(ID[ii])[0]) {
                     g(ID[ii]) += vtkGrid->GetPointData()->GetScalars("analytic_phi")->GetTuple(ID[ii])[0];
                 } else {
-                    g(ID[ii]) += (PI*V)*(vtkGrid->GetCellData()->GetScalars("G")->GetTuple(ic)[0]); // Why 3.0??
+                    g(ID[ii]) += (V/4.)*(vtkGrid->GetCellData()->GetScalars("G")->GetTuple(ic)[0]); // Why 3.0??
                 }
-                //if (std::abs(C(ii,3)-13.5) > 1e-5) {
-                //    g(ID[ii]) -= G[ii]*(V)*sigma_bar;
-                //}
-                //g(ID[ii]) +=  V/4*avg;
-                  //g(ID[ii]) +=  6.67 *(V/4.) * avg;
             }
-            //g(ID[0]) +=  (V/4.) * avg;
-        }
-        //g -= A*DB;
 
+        }
     }
 
-    void FEM4EllipticPDE::SolveOLD(const std::string& fname) {
+    void FEM4EllipticPDE::SolveOLD2(const std::string& fname) {
 
         Real r0[3];
         Real r1[3];