Adding linear mag code

examples/CMakeLists.txt

@@ -1,8 +1,8 @@
 add_executable( FEM4EllipticPDE_bhmag FEM4EllipticPDE_bhmag.cpp  )
 target_link_libraries(  FEM4EllipticPDE_bhmag  "lemmacore" "fem4ellipticpde")
 
-add_executable( FEM4EllipticPDE FEM4EllipticPDE.cpp  )
-target_link_libraries(  FEM4EllipticPDE  "lemmacore" "fem4ellipticpde")
+#add_executable( FEM4EllipticPDE FEM4EllipticPDE.cpp  )
+#target_link_libraries(  FEM4EllipticPDE  "lemmacore" "fem4ellipticpde")
 
 add_executable( merge merge.cpp  )
 target_link_libraries(  merge  "lemmacore" "fem4ellipticpde")
@@ -31,9 +31,8 @@ install(DIRECTORY "LinearMag"
 	PATTERN "CMakeLists.txt"  EXCLUDE
 )
 
-#INSTALL_TARGETS( "${CMAKE_INSTALL_PREFIX}/share/FEM4EllipticPDE/"
 INSTALL_TARGETS( "/share/FEM4EllipticPDE/"
-	FEM4EllipticPDE_bhmag FEM4EllipticPDE merge VTKDC VTKEdgeG VTKEdgeGsphere ResampleWithDataset
+	FEM4EllipticPDE_bhmag  merge VTKDC VTKEdgeG VTKEdgeGsphere ResampleWithDataset
 )
 
 INSTALL_TARGETS( "/share/FEM4EllipticPDE/grav"

examples/FEM4EllipticPDE_bhmag.cpp

@@ -93,7 +93,7 @@ int main(int argc, char**argv) {
         //std::cout << "usage: ./utFEMEllipticPDE_bhmag   <mesh.vtu>  <results.vtu>" << std::endl;
         exit(EXIT_SUCCESS);
     }
-
+/*
     int    nx = 80;
     //double dx =  .0021875 ; // 160
     double dx =  .004375 ; // 160
@@ -119,7 +119,7 @@ int main(int argc, char**argv) {
         rGrid->SetXCoordinates(xCoords);
         rGrid->SetYCoordinates(yCoords);
         rGrid->SetZCoordinates(zCoords);
-
+*/
     ////////////////////////////////////////////
     // Define Sigma/mu term
     // TODO
@@ -131,10 +131,10 @@ int main(int argc, char**argv) {
         GReader->SetFileName(argv[1]);
         GReader->Update();
 
-    vtkImplicitDataSet* implG = vtkImplicitDataSet::New();
-        implG->SetDataSet(GReader->GetOutput());
-        implG->SetOutValue(0.);
-        implG->SetOutGradient(0., 0., 0.);
+    //vtkImplicitDataSet* implG = vtkImplicitDataSet::New();
+    //    implG->SetDataSet(GReader->GetOutput());
+    //    implG->SetOutValue(0.);
+    //    implG->SetOutGradient(0., 0., 0.);
 
     ////////////////////////////////////////////
     // Define solution mesh
@@ -146,10 +146,10 @@ int main(int argc, char**argv) {
     ////////////////////////////////////////////
     // Solve
     FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
-        Solver->SetGFunction(implG);
+        //Solver->SetGFunction(implG);
         //Solver->SetGFunction(Magnet);
         //Solver->SetSigmaFunction(implSigma);
-        Solver->SetBoundaryStep(.05);
+        //Solver->SetBoundaryStep(.05);
         Solver->SetGrid(MeshReader->GetOutput());
         //Solver->SetupPotential();
         //Solver->SetGrid(rGrid);
@@ -158,7 +158,7 @@ int main(int argc, char**argv) {
     // Clean up
     Solver->Delete();
     GReader->Delete();
-    implG->Delete();
+    //implG->Delete();
 
     exit(EXIT_SUCCESS);
 }

examples/LinearMag/Sphere.cpp

@@ -18,12 +18,37 @@
  */
 
 #include "Lemma"
+#include "LinearMag.h"
 
 using namespace Lemma;
 
-int main() {
+int main( int argc, char** argv ) {
 
-    std::cout << "Test of LinearMag class\n";
 
+    std::cout << "###################################################\n";
+    std::cout << "#          Example of LinearMag class             #\n";
+    std::cout << "#                                                 #\n";
+    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\n";
+
+    if (argc < 3) {
+        std::cout << "\e[1museage\e[0m:\n";
+        std::cout << "./Sphere  in.vtu out.vtu\n\n" << std::endl;
+        exit(EXIT_SUCCESS);
+    }
+
+    LinearMag* Mag = LinearMag::New();
+        Mag->SetInducingMagField( 55234., 60, 14, NANOTESLA );
+        Mag->SetVTUGridFile( argv[1] );
+        Mag->CalculateRHS( "kappa" );
+        Mag->Solve( argv[2] );
+
+    // Print out YAML serialization of class
+    std::cout << *Mag << std::endl;
+
+    Mag->Delete();
+    exit(EXIT_SUCCESS);
 }
 

examples/VTKEdgeGsphere.cpp

@@ -81,8 +81,13 @@ int main(int argc, char**argv) {
     vtkDoubleArray* G = vtkDoubleArray::New();
         G->SetNumberOfComponents(1);
         G->SetNumberOfTuples(nn);
-        G->SetName("G");
-
+        //G->SetName("G");
+        if ( std::string(argv[3]) == "mag") {
+            G->SetName("kappa");
+        }
+        else {
+            G->SetName("G");
+        }
     vtkDoubleArray* phi = vtkDoubleArray::New();
         phi->SetNumberOfComponents(1);
         phi->SetNumberOfTuples(nn);
@@ -105,31 +110,25 @@ int main(int argc, char**argv) {
                 G->InsertTuple1( in, 0 );
             }
         }
-
-        else if ( std::string(argv[3]) == "magnet") {
+        if ( std::string(argv[3]) == "magnet") {
             if ( raddist > R - eps && raddist < R + eps ) {
                 // \rho = \nabla \cdot \mathbf{M}
-                //
-                if ( point[2] < -eps ) {
-                    G->InsertTuple1( in,   1 );
-                } else if (point[2] > eps) {
-                    G->InsertTuple1( in,  -1 );
-                } else {
-                    G->InsertTuple1( in, 0 );
-                }
-                //
-                // Above is inaccurate and unstable
-
                 // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
                 Eigen::Vector3d n;
                 n << point[0],point[1],point[2];
                 n.array() /= raddist;
-                //double sigma = n.dot(M);
                 G->InsertTuple1(in, n.dot(M) );
             } else {
                 G->InsertTuple1( in, 0 );
             }
         }
+        if ( std::string(argv[3]) == "mag") {
+            if ( raddist < R + eps) {
+                G->InsertTuple1( in, 1 );
+            } else {
+                G->InsertTuple1( in, 0 );
+            }
+        }
 
         // Insert analytic phi part
         /* magnetics problem p. 198 Jackson */

include/FEM4EllipticPDE.h

@@ -49,6 +49,8 @@
 #include "vtkDataSetSurfaceFilter.h"
 #include "vtkCellArray.h"
 #include "vtkCellData.h"
+#include <vtkXMLUnstructuredGridReader.h>
+#include <vtkUnstructuredGrid.h>
 
 #include <Eigen/IterativeLinearSolvers>
 #include <Eigen/SparseLU>
@@ -131,9 +133,23 @@ namespace Lemma {
              */
             void SetGFunction( Real (*pFcn3)(const Real&, const Real&, const Real&) );
 
+            /* Sets the vtkDataSet that defines the calculation grid.
+             */
+            //void SetGrid(vtkDataSet* Grid);
+
             /** Sets the vtkDataSet that defines the calculation grid.
              */
-            void SetGrid(vtkDataSet* Grid);
+            void SetGrid(vtkUnstructuredGrid* Grid);
+
+            /**
+             *  Read grid in from VTK file
+             */
+            void SetVTUGridFile( std::string& vtkGridFile );
+
+            /**
+             *  Read grid in from VTK file
+             */
+            void SetVTUGridFile( char* vtkGridFile );
 
             /** Sets up a DC problem with a survey
              *  @param[in] ij is the injection index
@@ -147,8 +163,8 @@ namespace Lemma {
             /** Performs solution */
             void Solve(const std::string& fname);
 
-            /** Performs solution */
-            void SolveOLD2(const std::string& fname);
+            /* Performs solution */
+            //void SolveOLD2(const std::string& fname);
 
             // ====================  ACCESS        =======================
 
@@ -289,7 +305,7 @@ namespace Lemma {
             vtkSmartPointer<vtkImplicitFunction>    vtkG;
 
             /** Any type of vtkDataSet may be used as a calculation Grid */
-            vtkSmartPointer<vtkDataSet>             vtkGrid;
+            vtkSmartPointer<vtkUnstructuredGrid>             vtkGrid;
 
             /** Function pointer to function describing g */
             Real (*gFcn3)(const Real&, const Real&, const Real&);

include/LinearMag.h

@@ -20,7 +20,8 @@
 #ifndef  LINEARMAG_INC
 #define  LINEARMAG_INC
 
-#include	"FEM4EllipticPDE.h"
+#include  <vtkCellIterator.h>
+#include  "FEM4EllipticPDE.h"
 
 namespace Lemma {
 
@@ -78,6 +79,28 @@ namespace Lemma {
 
         // ====================  ACCESS        =======================
 
+        /**
+         *  @param[in] B0 is the incident magnetic (-induction) field
+         *  @param[in] U is a MAGUNITS enum specifying the units of measurement
+         */
+        void SetInducingMagFieldVector( const Vector3r& B0, const MAGUNITS& U );
+
+        /**
+         *  @param[in] intensity is the incident magnetic (-induction) field intensity
+         *  @param[in] inc is the incident magnetic (-induction) field inclination in degrees
+         *  @param[in] dec is the incident magnetic (-induction) field declination in degrees
+         *  @param[in] U is a MAGUNITS enum specifying the units of intensity
+         */
+        void SetInducingMagField( const Real& intensity, const Real& inclination,
+                const Real& declination, const MAGUNITS& U );
+
+
+        /**
+         *  Calculates the right hand side of the equation
+         *   \f$ \nabla \cdot \kappa \mathbf{H}_0 \f$
+         */
+        void CalculateRHS( const std::string& susName );
+
         // ====================  INQUIRY       =======================
 
         protected:
@@ -105,6 +128,14 @@ namespace Lemma {
 
         // ====================  DATA MEMBERS  =========================
 
+        Vector3r    B0;
+
+
+        /**
+         *   Used internally to scale mag units into Tesla
+         */
+        void ScaleB0 ( const MAGUNITS& U);
+
     }; // -----  end of class  LinearMag  -----
 
     /*! @} */ // End of group

src/FEM4EllipticPDE.cpp

@@ -132,10 +132,45 @@ namespace Lemma {
         gFcn3 = gFcn;
     }
 
-    void FEM4EllipticPDE::SetGrid(vtkDataSet* grid) {
+    //void FEM4EllipticPDE::SetGrid(vtkDataSet* grid) {
+    //    vtkGrid = grid;
+    //}
+
+    void FEM4EllipticPDE::SetGrid(vtkUnstructuredGrid* grid) {
         vtkGrid = grid;
     }
 
+    //--------------------------------------------------------------------------------------
+    //       Class:  FEM4EllipticPDE
+    //      Method:  SetVTUGridFile
+    //--------------------------------------------------------------------------------------
+    void FEM4EllipticPDE::SetVTUGridFile ( std::string& fname  ) {
+        vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
+        MeshReader->SetFileName(fname.c_str());
+        MeshReader->Update();
+        //vtkGrid->DeepCopy( MeshReader->GetOutput() );
+        vtkGrid->ShallowCopy( MeshReader->GetOutput() );
+        MeshReader->Delete();
+        return ;
+    }		// -----  end of method FEM4EllipticPDE::SetVTKGridFile  -----
+
+    //--------------------------------------------------------------------------------------
+    //       Class:  FEM4EllipticPDE
+    //      Method:  SetVTUGridFile
+    //--------------------------------------------------------------------------------------
+    void FEM4EllipticPDE::SetVTUGridFile ( char* fname  ) {
+        std::cout  << "Loading VTK .vtu file " << fname << std::endl;
+        vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
+        MeshReader->SetFileName(fname);
+        MeshReader->Update();
+        //vtkGrid->DeepCopy( MeshReader->GetOutput() );
+        vtkGrid = MeshReader->GetOutput();
+        //vtkGrid->ShallowCopy( MeshReader->GetOutput() );
+        MeshReader->Delete();
+        return ;
+    }		// -----  end of method FEM4EllipticPDE::SetVTKGridFile  -----
+
+
     vtkSmartPointer<vtkIdList> FEM4EllipticPDE::GetConnectedPoints(const int& id0) {
         vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
         vtkSmartPointer<vtkIdList> cellList = vtkSmartPointer<vtkIdList>::New();
@@ -549,10 +584,7 @@ namespace Lemma {
             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
+            Real i4pi = .5;
             /* surfaces of tetrahedra */
             // Face 0, ID 0,1,2
             Eigen::Matrix<Real, 3, 3> CC = Eigen::Matrix<Real, 3, 3>::Ones() ;
@@ -602,6 +634,7 @@ namespace Lemma {
         }
     }
 
+    #if 0
     void FEM4EllipticPDE::SolveOLD2(const std::string& fname) {
 
         Real r0[3];
@@ -888,6 +921,7 @@ namespace Lemma {
         uArray->Delete();
 
     }
+    #endif
 
     // Uses simpon's rule to perform a definite integral of a
     // function (passed as a pointer). The integration occurs from

src/LinearMag.cpp

@@ -104,6 +104,7 @@ YAML::Node  LinearMag::Serialize (  ) const {
     YAML::Node node = FEM4EllipticPDE::Serialize();;
     node.SetTag( this->Name );
     // FILL IN CLASS SPECIFICS HERE
+    node["B0"] = B0;
     return node;
 }		// -----  end of method LinearMag::Serialize  -----
 
@@ -120,4 +121,111 @@ LinearMag* LinearMag::DeSerialize ( const YAML::Node& node  ) {
 }		// -----  end of method LinearMag::DeSerialize  -----
 #endif
 
+
+
+//--------------------------------------------------------------------------------------
+//       Class:  LinearMag
+//      Method:  SetInducingMagField
+//--------------------------------------------------------------------------------------
+void LinearMag::SetInducingMagField ( const Real& intensity, const Real& inc,
+                const Real& dec, const MAGUNITS& U ) {
+
+    B0(0) = intensity * std::cos(inc) * std::cos(dec); // northing
+    B0(1) = intensity * std::cos(inc) * std::sin(dec); // easting
+    B0(2) = intensity * std::sin(inc);                 // z
+    ScaleB0(U);
+    return ;
+}		// -----  end of method LinearMag::SetInducingMagField  -----
+
+
+//--------------------------------------------------------------------------------------
+//       Class:  LinearMag
+//      Method:  SetInducingMagFieldVector
+//--------------------------------------------------------------------------------------
+void LinearMag::SetInducingMagFieldVector ( const Vector3r& BB0, const MAGUNITS& U  ) {
+    B0 = BB0;
+    ScaleB0(U);
+    return ;
+}		// -----  end of method LinearMag::SetInducingMagFieldVector  -----
+
+
+//--------------------------------------------------------------------------------------
+//       Class:  LinearMag
+//      Method:  ScaleB0
+//--------------------------------------------------------------------------------------
+void LinearMag::ScaleB0 ( const MAGUNITS& U ) {
+    switch ( U ) {
+        case TESLA:
+            break;
+        case NANOTESLA:
+            B0 *= 1e-9;
+            break;
+        case GAUSS:
+            B0 *= 1e-4;
+            break;
+    }
+    return ;
+}		// -----  end of method LinearMag::ScaleB0  -----
+
+
+//--------------------------------------------------------------------------------------
+//       Class:  LinearMag
+//      Method:  CalculateRHS
+//--------------------------------------------------------------------------------------
+void LinearMag::CalculateRHS ( const std::string& susName ) {
+    vtkDoubleArray* G = vtkDoubleArray::New();
+    G->SetNumberOfComponents(1);
+    G->SetNumberOfTuples( vtkGrid->GetNumberOfPoints() );
+    G->SetName("G");
+
+    // Iterate over all the points or all of the cells?
+    for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {
+
+        Real cellSus = vtkGrid->GetCellData()->GetScalars(susName.c_str())->GetTuple(ic)[0];
+
+        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);
+            C(ii, 0) = 1;
+            C(ii, 1) =  pts[0];
+            C(ii, 2) =  pts[1];
+            C(ii, 3) =  pts[2];
+        }
+
+        /* 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);
+
+        /* 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]
+        */
+        // 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 ;
+        }
+        */
+    }
+
+    vtkGrid->GetPointData()->AddArray( G );
+    return ;
+}		// -----  end of method LinearMag::CalculateRHS  -----
+
+
 }		// -----  end of Lemma  name  -----