Appears to be working with Dirichlet (non-homogeneous) conditions

examples/CMakeLists.txt

@@ -16,12 +16,15 @@ 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( 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 VTKEdgeGsphere ResampleWithDataset
+	FEM4EllipticPDE_bhmag FEM4EllipticPDE merge VTKDC VTKEdgeG VTKGsphere VTKEdgeGsphere ResampleWithDataset
 )
 
 install( DIRECTORY "borehole"

examples/VTKGsphere.cpp

@@ -0,0 +1,164 @@
+/* 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/.
+ */
+
+/**
+ * @file
+ * @date      08/07/2014 04:16:25 PM
+ * @version   $Id$
+ * @author    Trevor Irons (ti)
+ * @email     Trevor.Irons@xri-geo.com
+ * @copyright Copyright (c) 2014, XRI Geophysics, LLC
+ * @copyright Copyright (c) 2014, Trevor Irons
+ */
+
+#include <vtkUnstructuredGrid.h>
+
+#include <vtkUnstructuredGridReader.h>
+#include <vtkUnstructuredGridWriter.h>
+
+#include <vtkGenericDataObjectReader.h>
+
+#include <vtkXMLUnstructuredGridReader.h>
+#include <vtkXMLUnstructuredGridWriter.h>
+#include <vtkDoubleArray.h>
+#include <vtkPointData.h>
+#include <vtkCellData.h>
+#include <vtkCell.h>
+#include <vtkCellCenters.h>
+#include <cmath>
+#include <Eigen/Core>
+
+const double PI = 4.0*atan(1.0);
+
+int main(int argc, char**argv) {
+
+    std::cout << "Mesh processing routine\n";
+    if (argc<4) {
+        std::cout << "usage:\n" << "./utVTKEdge  inMesh.vtu  <radius> <problemType>  outG.vtu" << std::endl;
+        exit(EXIT_SUCCESS);
+    }
+
+    vtkUnstructuredGrid* uGrid = vtkUnstructuredGrid::New();
+
+    std::string fn = argv[1];
+    if(fn.substr(fn.find_last_of(".") + 1) == "vtk") {
+        vtkGenericDataObjectReader* Reader = vtkGenericDataObjectReader::New();
+            Reader->SetFileName(argv[1]);
+            Reader->Update();
+        if(Reader->IsFileUnstructuredGrid()) {
+            std::cout << "Found Unstructured grid legacy file" << std::endl;
+            uGrid = Reader->GetUnstructuredGridOutput();
+        } else {
+            std::cerr << "Unknown legacy format";
+            exit(EXIT_FAILURE);
+        }
+
+    } else {
+
+        vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
+            std::cout << "Reading" << argv[1] << std::endl;
+            Reader->SetFileName(argv[1]);
+            Reader->Update();
+            uGrid = Reader->GetOutput();
+    }
+
+
+        int nc = uGrid->GetNumberOfCells()  ;
+        int nn = uGrid->GetNumberOfPoints()  ;
+
+    std::cout << "Processing grid with nodes=" << nn << "\t cells=" << nc << "\n";
+
+    // Input magnet size TODO get from file or .geo file even?
+    double R =  atof(argv[2]); //.25;          // radius of magnet
+    double eps = 1e-4;       // epsilon for on the point
+    // Pol = double[3];
+
+    // Declare new array for the data
+    vtkDoubleArray* G = vtkDoubleArray::New();
+        G->SetNumberOfComponents(1);
+        G->SetNumberOfTuples(nc);
+        G->SetName("G");
+
+    vtkDoubleArray* phi = vtkDoubleArray::New();
+        phi->SetNumberOfComponents(1);
+        phi->SetNumberOfTuples(nn);
+        phi->SetName("analytic_phi");
+
+    // Set point data
+    double point[3];
+    Eigen::Vector3d M; M << 0,0,1;
+    for (int in=0; in<nn; ++in) {
+        uGrid->GetPoint(in, &point[0]);
+        //std::cout << point[0] << "\t" <<point[1] << "\t" << point[2] << std::endl;
+        double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
+
+        //double rho = sqrt( point[1]*point[1] + point[2]*point[2] );
+        // Calculate RHS
+
+        // Insert analytic phi part
+        /* magnetics problem p. 198 Jackson */
+        if (std::string(argv[3]) == "magnet") {
+            if (raddist < R) {
+                phi->InsertTuple1( in, (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, (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)
+            if (raddist < R) {
+                phi->InsertTuple1( in, mass * (( 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 );
+            }
+        }
+    }
+
+    vtkCellCenters* cellCentersFilter = vtkCellCenters::New();
+    cellCentersFilter->SetInputData( uGrid );
+    cellCentersFilter->VertexCellsOn();
+    cellCentersFilter->Update();
+
+    // Set Cell Data
+    for (int ic=0; ic<nc; ++ic) {
+        //uGrid->GetCell(ic)->GetParametricCenter( &point[0] );
+        cellCentersFilter->GetOutput()->GetPoint(ic, &point[0]);
+        double raddist = sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] );
+        //std::cout << "point " << point[0] << "\t" << point[1] << "\t" << point[2] << std::endl;
+        if ( std::string(argv[3]) == "gravity") {
+            if ( raddist < R + eps) {
+                G->InsertTuple1( ic, 1 );
+            } else {
+                G->InsertTuple1( ic, 0 );
+            }
+        }
+
+    }
+
+    // Add new data
+    uGrid->GetPointData()->AddArray( phi );
+    uGrid->GetCellData()->AddArray( G );
+
+    // Write out new file
+    vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
+        //Writer->SetInputConnection(Reader->GetOutputPort());
+        Writer->SetInputData( uGrid );
+        Writer->SetFileName( argv[4] );
+        Writer->Write();
+
+    //Reader->Delete();
+
+
+}

examples/borehole/sphere.geo

@@ -7,11 +7,11 @@
  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
  */
 
-radius = 2.25;     // Radius of the damn thing
-lc = radius/4;     //  0.25;   // Target element size
+radius = 4.25;     // Radius of the damn thing
+lc = radius/5;     //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 3*radius; // The down side of potential
+Box = 5*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;
@@ -19,7 +19,7 @@ Y1 =  Box;
 Z0 = -Box;
 Z1 =  Box;
 
-cellSize=radius/32; ///10;
+cellSize=radius/10; ///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 
-../VTKEdgeGsphere MergedSphere.vtu   2.25  gravity  sphereGrav.vtu
+../VTKGsphere MergedSphere.vtu   4.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 = 2.25;     // Radius of the damn thing
+radius = 4.25;     // Radius of the damn thing
 
 lc = radius;   //  0.25;   // Target element size
 
 // Total Solution Space
-Box = 3*radius; // The down side of potential
+Box = 5*radius; // The down side of potential
 X0 = -Box;
 X1 =  Box;
 Y0 = -Box;

src/FEM4EllipticPDE.cpp

@@ -436,7 +436,7 @@ 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]) += (3.0*V)*(vtkGrid->GetPointData()->GetScalars("G")->GetTuple(ID[ii])[0]); // Why 3.0??
+                    g(ID[ii]) += (PI*V)*(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;