Exmple for Andy and Boxin

examples/borehole/sphere.geo

 lc = radius;   //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 10*radius; // The down side of potential
+Box = 30*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;
 ///////////////////////////////////////////////
 // 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("(2.25 - F1)^2 + %g", cellSize*10 );  // WORKS

examples/borehole/sphere.sh

-
+#!/usr/bin/env bash
 gmsh  -3 -format vtk -o sphere.vtk  sphere.geo 
 gmsh  -2 -format stl -o sphereBox.stl  sphereBox.geo 
+paraview --state=boundarySphere.pvsm 
 
+#IMPORTANT! 
+#   Select ResampleWithDataset1
+#	-->File -> SaveData -> MergedSphere.vtu 
 
+# Set up problem 
 
-#gmsh -3 -format vtk -o sphere.vtk  sphere.geo
-#../utVTKEdgeGsphere sphere.vtk .25 magnet  sphereG.vtu
-#../utFEM4EllipticPDE_bhmag  sphereG.vtu  sphereG.vtu   sphereOut.vtu 
 
 # Gravity problem
-../utVTKEdgeGsphere sphere.vtk .25 gravity  sphereGrav.vtu
-../utFEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
+#../VTKEdgeGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
+#../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
+
 #paraview --state=sliceSphere.pvsm
 
 

examples/borehole/sphere2.sh

+#!/usr/bin/env bash
+#gmsh  -3 -format vtk -o sphere.vtk  sphere.geo 
+#gmsh  -2 -format stl -o sphereBox.stl  sphereBox.geo 
+#paraview --state=boundarySphere.pvsm 
+
+#IMPORTANT! 
+#   Select ResampleWithDataset1
+#	-->File -> SaveData -> MergedSphere.vtu 
+
+# Set up problem 
+# Gravity problem
+../VTKEdgeGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
+../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
+
+#paraview --state=sliceSphere.pvsm
+
+
+# --data=sphereOut.vtu
+
+# even finer meshing
+# radius = 4, u = \pm .9              analytic = \pm 1.33
+
+# Finer meshing
+# radius = 5, u = \pm 1.65            analytic = \pm 1.67
+# radius = 4, u = \pm 1.09            analytic = \pm 1.33
+
+# radius = 5, u = \pm 2.5             analytic = \pm 1.67
+# radius = 4, u = \pm 1.66            analytic = \pm 1.33
+# radius = 3, u= \pm .926 -.985       analytic = \pm 1
+# radius = 2, u = \pm.42              analytic = \pm .667
+# for radius=1,  u = \pm .111         analytic = \pm .333 
+# for radius=.5, u = \pm .0477 -.0462 analytic = \pm .167

examples/borehole/sphereBox.geo

 lc = radius;   //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 10*radius; // The down side of potential
+Box = 30*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;

src/FEM4EllipticPDE.cpp

         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];
+//             for (int ii=0; ii<4; ++ii) {
+//                 double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
+//                 C(ii, 0) = 1;
+//                 C(ii, 1) =  pts[0];
+//                 C(ii, 2) =  pts[1];
+//                 C(ii, 3) =  pts[2];
 //             }
 
 //             Real V = (1./6.) * C.determinant();          // volume of tetrahedra
         //Eigen::BiCGSTAB<Eigen::SparseMatrix<Real> > cg(A);
         cg.setMaxIterations(3000);
 
+        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;
 
         // Build stiffness matrix (A)
         std::cout << "Building Stiffness Matrix (A) " << std::endl;
         std::cout << "\tMesh has " << vtkGrid->GetNumberOfCells()  << " cells "  << std::endl;
+        std::cout << "\tMesh has " << vtkGrid->GetNumberOfPoints()  << " points "  << std::endl;
 
   	    //Eigen::SparseMatrix<Real>
         A.resize(vtkGrid->GetNumberOfPoints(), vtkGrid->GetNumberOfPoints());
 
             // construct coordinate matrix C
             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] ;
+            for (int ii=0; ii<4; ++ii) {
+                double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
+                C(ii, 0) = 1;
+                C(ii, 1) =  pts[0] ;
+                C(ii, 2) =  pts[1] ;
+                C(ii, 3) =  pts[2] ;
             }
 
             Eigen::Matrix<Real, 4, 4> Phi = C.inverse(); // nabla \phi
                 ID[1] = Ids->GetId(1);
                 ID[2] = Ids->GetId(2);
                 ID[3] = Ids->GetId(3);
-            Real sum(0), sigma_bar(0);
+
+            Real sigma_bar(0);
             if (GCell) {
                 sigma_bar = vtkGrid->GetCellData()->GetScalars("G")->GetTuple1(ic);
             } else {
                 sigma_bar += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3]);
                 sigma_bar /= 4.;
             }
+            sigma_bar = 1.;
 
-            for (int ip=0; ip<4; ++ip) {
-                for (int ip2=0; ip2<4; ++ip2) {
-                    if (ip2 == ip) {
+            for (int ii=0; ii<4; ++ii) {
+                for (int jj=0; jj<4; ++jj) {
+                    if (jj == ii) {
                         // I apply boundary to Stiffness matrix, it's common to take the other approach and apply to the load vector and then
                         //   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];
+                        //Real bdry = V*(1./(BndryH*BndryH))*BndrySigma*bndryCnt( ID[ii] ); // + sum;
+                        Real bb = vtkGrid->GetPointData()->GetScalars("vtkValidPointMask")->GetTuple(ID[ii])[0];
                         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 ) );
+                        coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], bdry +  Phi.col(ii).tail<3>().dot(Phi.col(jj).tail<3>() ) * V * sigma_bar ) );
                     } else {
-                        coeffs.push_back( Eigen::Triplet<Real> ( ID[ip], ID[ip2],         Phi.col(ip).tail<3>().dot(Phi.col(ip2).tail<3>() ) * V * sigma_bar ) );
+                        coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj],         Phi.col(ii).tail<3>().dot(Phi.col(jj).tail<3>() ) * V * sigma_bar ) );
                     }
                     // Stiffness matrix no longer contains boundary conditions...
-                    //coeffs.push_back( Eigen::Triplet<Real> ( ID[ip], ID[ip2],         Phi.col(ip).tail<3>().dot(Phi.col(ip2).tail<3>() ) * V * sigma_bar ) );
+                    //coeffs.push_back( Eigen::Triplet<Real> ( ID[ii], ID[jj], Phi.col(ii).tail<3>().dot(Phi.col(jj).tail<3>() ) * V * sigma_bar ) );
                 }
             }
             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;
     }
 
     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;
+        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 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];
-                }
+            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] ;
+                C(ii, 0) = 1;
+                C(ii, 1) =  pts[0];
+                C(ii, 2) =  pts[1];
+                C(ii, 3) =  pts[2];
+            }
 
-                Real V = (1./6.) * C.determinant();          // volume of tetrahedra
-                //Real V = C.determinant();          // volume of tetrahedra
+            Real V = (1./6.) * 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);
 
-                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 )
+            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];
             }
 
+            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.;
+                //g(ID[ii]) +=  (V/4.) * ( vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0] )  ;
+                g(ID[ii]) +=  (V/4.) * avg;
+            }
         }
 
     }
         // Expensive search for whether or not point is on boundary. O(n) cost.
         VectorXi bndryCnt = VectorXi::Zero(vtkGrid->GetNumberOfPoints());
         for (int isp=0; isp < Surface->GetOutput()->GetNumberOfPoints(); ++isp) {
-            //double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ip); //[ipc] ;
+            //double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
             // x \in -14.5 to 14.5
             // y \in 0 to 30
             bndryCnt(BdryIds->GetTuple1(isp)) += 1;
         // Build stiffness matrix (A)
         std::cout << "Building Stiffness Matrix (A) " << std::endl;
         std::cout << "\tMesh has " << vtkGrid->GetNumberOfCells()  << " cells "  << std::endl;
+        std::cout << "\tMesh has " << vtkGrid->GetNumberOfPoints()  << " points "  << std::endl;
 
 
   	    Eigen::SparseMatrix<Real> A(vtkGrid->GetNumberOfPoints(), vtkGrid->GetNumberOfPoints());
 
             // construct coordinate matrix C
             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] ;
+            for (int ii=0; ii<4; ++ii) {
+                double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
+                C(ii, 0) = 1;
+                C(ii, 1) =  pts[0] ;
+                C(ii, 2) =  pts[1] ;
+                C(ii, 3) =  pts[2] ;
             }
 
-            Eigen::Matrix<Real, 4, 4> Phi = C.inverse(); // nabla \phi
+            Eigen::Matrix<Real, 4, 4> GradPhi = C.inverse(); // nabla \phi
             Real V = (1./6.) * C.determinant();          // volume of tetrahedra
 
             vtkIdList* Ids = vtkGrid->GetCell(ic)->GetPointIds();
             Real sum(0);
             Real sigma_bar = vtkGrid->GetCellData()->GetScalars()->GetTuple1(ic);
 
-            for (int ip=0; ip<4; ++ip) {
-                for (int ip2=0; ip2<4; ++ip2) {
-                    if (ip2 == ip) {
+            for (int ii=0; ii<4; ++ii) {
+                for (int jj=0; jj<4; ++jj) {
+                    if (jj == ii) {
                         // I apply boundary to Stiffness matrix, it's common to take the other approach and apply to the load vector and then
                         //   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];
+                        //Real bdry = V*(1./(BndryH*BndryH))*BndrySigma*bndryCnt( ID[ii] ); // + sum;
+                        Real bb = vtkGrid->GetPointData()->GetScalars("vtkValidPointMask")->GetTuple(ID[ii])[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 ) );
+                        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 ) );
                     } else {
-                        coeffs.push_back( Eigen::Triplet<Real> ( ID[ip], ID[ip2],         Phi.col(ip).tail<3>().dot(Phi.col(ip2).tail<3>() ) * V * sigma_bar ) );
+                        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[ip], ID[ip2],         Phi.col(ip).tail<3>().dot(Phi.col(ip2).tail<3>() ) * V * sigma_bar ) );
+                    //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 ;
             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];
+                for (int ii=0; ii<4; ++ii) {
+                    double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
+                    C(ii, 0) = 1;
+                    C(ii, 1) =  pts[0];
+                    C(ii, 2) =  pts[1];
+                    C(ii, 3) =  pts[2];
                 }
 
                 Real V = (1./6.) * C.determinant();          // volume of tetrahedra
                 /*
                 if (!pe) {
                     if (std::abs(pt[0]) < .2 && std::abs(pt[1]-15) < .2 && pt[2] < 3.25 ) {
-                        sum = 1; //vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ip])[0];
+                        sum = 1; //vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0];
                         pe = true;
                     }
                 }*/
 /*
                 if (!ne) {
                     if (std::abs(pt[0]+1.) < .2 && std::abs(pt[1]-15) < .2 && pt[2] < 3.25 ) {
-                        sum = -1; //vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ip])[0];
+                        sum = -1; //vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0];
                         std::cout << "Negative Electroce\n";
                         ne = true;
                     }
                 }
 */
-                //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 )
+                //for (int ii=0; ii<4; ++ii) {
+                    //g(ID[ii]) +=  (V/4.) *  ( vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0] )  ;
+                    //if ( std::abs(vtkGrid->GetPointData()->GetScalars()->GetTuple(ID[ii])[0]) > 1e-3 )
                 //}
                 // TODO check Load Vector...
                 g(ID[0]) = sum; //(V/4.) * sum; // Why 3, Why!!!?
 
                     for(vtkIdType i = 0; i < connectedVertices->GetNumberOfIds(); ++i) {
 
-                        int ip = connectedVertices->GetId(i);
-                        vtkGrid->GetPoint(ip, r1);
-                        double g1 = vtkGrid->GetPointData()->GetScalars()->GetTuple(ip)[0] ;
+                        int ii = connectedVertices->GetId(i);
+                        vtkGrid->GetPoint(ii, r1);
+                        double g1 = vtkGrid->GetPointData()->GetScalars()->GetTuple(ii)[0] ;
                         //g(ic) += g0*dist(r0,r1); //CompositeSimpsons2(g0, r0, r1, 8);
                         if ( std::abs(g1) > 1e-3 ) {
                             g(ic) += CompositeSimpsons2(g1, g0, r1, r0, 1000);
                 //g(ic) += vtkGrid->GetPointData()->GetScalars()->GetTuple1(ic);// FunctionValue(r0[0], r0[1], r0[2]) ;
                 /*
                 for(vtkIdType i = 0; i < connectedVertices->GetNumberOfIds(); ++i) {
-                    int ip = connectedVertices->GetId(i);
-                    vtkGrid->GetPoint(ip, r1);
+                    int ii = connectedVertices->GetId(i);
+                    vtkGrid->GetPoint(ii, r1);
                     g(ic) += CompositeSimpsons2(vtkG, r0, r1, 8); // vtkG->FunctionValue(r0[0], r0[1], r0[2]) ;
                 }
                 */
                 //#pragma omp parallel for reduction(+:sum)
                 //#endif
                 for(vtkIdType i = 0; i < connectedVertices->GetNumberOfIds(); ++i) {
-                    int ip = connectedVertices->GetId(i);
-                    vtkGrid->GetPoint(ip, r1);
+                    int ii = connectedVertices->GetId(i);
+                    vtkGrid->GetPoint(ii, r1);
                     g(ic) += CompositeSimpsons2(gFcn3, r0, r1, 8); // vtkG->FunctionValue(r0[0], r0[1], r0[2]) ;
                     //sum += CompositeSimpsons2(gFcn3, r0, r1, 8); // vtkG->FunctionValue(r0[0], r0[1], r0[2]) ;
                 }