Surface sources now working, but need to add switch mechanisms.

examples/CMakeLists.txt

@@ -16,15 +16,19 @@ target_link_libraries( VTKEdgeG  "lemmacore" "fem4ellipticpde")
 add_executable( VTKEdgeGsphere VTKEdgeGsphere.cpp  )
 target_link_libraries( VTKEdgeGsphere  "lemmacore" "fem4ellipticpde")
 
-add_executable( VTKGsphere VTKGsphere.cpp  )
-target_link_libraries( VTKGsphere  "lemmacore" "fem4ellipticpde")
+add_executable( setupGravSphere setupGravSphere.cpp  )
+target_link_libraries( setupGravSphere  "lemmacore" "fem4ellipticpde")
 
 add_executable( ResampleWithDataset ResampleWithDataset.cpp  )
 target_link_libraries( ResampleWithDataset  "lemmacore" "fem4ellipticpde")
 
 #INSTALL_TARGETS( "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
 INSTALL_TARGETS( "/share/FEM4EllipticPDE/"
-	FEM4EllipticPDE_bhmag FEM4EllipticPDE merge VTKDC VTKEdgeG VTKGsphere VTKEdgeGsphere ResampleWithDataset
+	FEM4EllipticPDE_bhmag FEM4EllipticPDE merge VTKDC VTKEdgeG VTKEdgeGsphere ResampleWithDataset
+)
+
+INSTALL_TARGETS( "/share/FEM4EllipticPDE/grav"
+	setupGravSphere 
 )
 
 install( DIRECTORY "borehole"
@@ -32,6 +36,16 @@ install( DIRECTORY "borehole"
 	PATTERN .svn EXCLUDE
 )
 
+install( DIRECTORY "grav"
+	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
+	PATTERN .svn EXCLUDE
+)
+
+install( DIRECTORY "magnet"
+	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
+	PATTERN .svn EXCLUDE
+)
+
 install( DIRECTORY "DC"
 	DESTINATION "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
 	PATTERN .svn EXCLUDE

examples/borehole/sphere.sh → examples/grav/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   2.25  gravity  sphereGrav.vtu
+./setupGravSphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
 ../FEM4EllipticPDE_bhmag  sphereGrav.vtu  sphereGrav.vtu   sphereOutGrav.vtu 
 
 #paraview --data=sphereOutGrav.vtu 

include/FEM4EllipticPDE.h

@@ -228,6 +228,22 @@ namespace Lemma {
              */
             Real dist(const Vector3r& r0, const Vector3r& r1);
 
+            /** Sets up the potential problem from the VTK file
+             */
+            void SetupLineSourcePotential();
+
+            /** Sets up the potential problem from the VTK file
+             */
+            void SetupSurfaceSourcePotential();
+
+            /** Sets up the potential problem from the VTK file
+             */
+            void SetupPointSourcePotential();
+
+            /** Sets up the potential problem from the VTK file
+             */
+            void SetupVolumeSourcePotential();
+
             // ====================  DATA MEMBERS  =========================
 
             /** The effective step at the boundary condition. Defaults to 1

src/FEM4EllipticPDE.cpp

@@ -201,7 +201,9 @@ namespace Lemma {
 
     void FEM4EllipticPDE::Solve( const std::string& resfile ) {
         ConstructAMatrix();
-        SetupPotential();
+        //SetupLineSourcePotential();
+        SetupSurfaceSourcePotential();
+        //SetupPotential();
         //ConstructLoadVector();
 
         std::cout << "\nSolving" << std::endl;
@@ -389,10 +391,164 @@ namespace Lemma {
                 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]) += (V/4.)*(vtkGrid->GetCellData()->GetScalars("G")->GetTuple(ic)[0]); // Why 3.0??
+                    g(ID[ii]) += (V/4.)*(vtkGrid->GetCellData()->GetScalars("G")->GetTuple(ic)[0]); // Cell data
+                    /* for point data, determine if it is a point or line source */
+                    //g(ID[ii]) = (vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0]); // Point source
+                    //g(ID[ii]) = (vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0]); // Point source
                 }
             }
 
+
+
+        }
+    }
+
+    void FEM4EllipticPDE::SetupLineSourcePotential() {
+
+        std::cerr << " FEM4EllipticPDE::SetupLineSourcePotential is not completed\n";
+        exit(1);
+
+        ////////////////////////////////////////////////////////////
+	    // 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] ;
+                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::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);
+                ID[1] = Ids->GetId(1);
+                ID[2] = Ids->GetId(2);
+                ID[3] = Ids->GetId(3);
+
+            /* Line source between nodes */
+            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];
+                } else {
+                    Real nv1 = vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0];
+                    for (int jj=ii+1; jj<4; ++jj) {
+                        Real nv2 = vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[jj])[0];
+                        if ( std::abs(nv1) > 1e-12 && std::abs(nv2) > 1e-12) {
+                            Real length = ( C.row(ii).tail<3>()-C.row(jj).tail<3>()  ).norm();
+
+                            g(ID[ii]) += ((nv1+nv2)/2.)/(length);//*(V/4.);// * (vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0]); // Point source
+                        }
+                    }
+                }
+            }
+        }
+
+    }
+
+    void FEM4EllipticPDE::SetupSurfaceSourcePotential() {
+
+        ////////////////////////////////////////////////////////////
+	    // 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] ;
+                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::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);
+                ID[1] = Ids->GetId(1);
+                ID[2] = Ids->GetId(2);
+                ID[3] = Ids->GetId(3);
+
+            // 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];
+                }
+            }
+
+            /* the 4 faces of the tetrahedra
+            ID[0] ID[1] ID[2]
+            ID[0] ID[1] ID[3]
+            ID[0] ID[2] ID[3]
+            ID[1] ID[2] ID[3]
+            */
+            //Real i4pi = .99738684646226885559*.961*PI/6.;//(4.*PI);// * V/4.;
+            //Real i4pi = .99738684646226885559*.961*PI/6.;//(4.*PI);// * V/4.;
+            //Real i4pi =  0.9584887594502403*PI/6.;//(4.*PI);// * V/4.;
+            Real i4pi = .5;// 0.1597132348845018*PI; // very nearly .5
+            /* surfaces of tetrahedra */
+            // Face 0, ID 0,1,2
+            Eigen::Matrix<Real, 3, 3> CC = Eigen::Matrix<Real, 3, 3>::Ones() ;
+            if (    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])) > 1e-12    ) {
+                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>();
+                Real TA = .5*CC.col(1).cross(CC.col(2)).norm();
+                g(ID[0]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0])*TA/3. * i4pi ;
+                g(ID[1]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1])*TA/3. * i4pi ;
+                g(ID[2]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])*TA/3. * i4pi ;
+            }
+            // Face 1, ID 0,1,3
+            if (    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])) > 1e-12    ) {
+                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>();
+                Real TA = .5*CC.col(1).cross(CC.col(2)).norm();
+                g(ID[0]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0])*TA/3. * i4pi ;
+                g(ID[1]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1])*TA/3. * i4pi ;
+                g(ID[3]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])*TA/3. * i4pi ;
+            }
+            // Face 2, ID 0,2,3
+            if (    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])) > 1e-12    ) {
+                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>();
+                Real TA = .5*CC.col(1).cross(CC.col(2)).norm();
+                g(ID[0]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[0])*TA/3. * i4pi ;
+                g(ID[2]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])*TA/3. * i4pi ;
+                g(ID[3]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])*TA/3. * i4pi ;
+            }
+            // Face 3, ID 1,2,3
+            if (    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])) > 1e-12 &&
+                    std::abs(vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])) > 1e-12    ) {
+                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>();
+                Real TA = .5*CC.col(1).cross(CC.col(2)).norm();
+                g(ID[1]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[1])*TA/3. * i4pi ;
+                g(ID[2]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[2])*TA/3. * i4pi ;
+                g(ID[3]) += vtkGrid->GetPointData()->GetScalars("G")->GetTuple1(ID[3])*TA/3. * i4pi ;
+            }
         }
     }