Added example files but no CMakeLists.txt yet. Testing branching.

examples/utFEM4EllipticPDE.cpp

@@ -0,0 +1,493 @@
+// ===========================================================================
+//
+//       Filename:  utFEM4EllipticPDE.cpp
+//
+//        Created:  08/16/12 19:49:10
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     08/16/12
+  @version   0.0
+ **/
+
+#include "FEM4EllipticPDE.h"
+#include "vtkSphere.h"
+#include "vtkRectilinearGrid.h"
+#include "vtkFloatArray.h"
+#include <vtkIdList.h>
+#include <vtkIdTypeArray.h>
+#include <vtkCell.h>
+#include <vtkTriangleFilter.h>
+#include <vtkDataSetSurfaceFilter.h>
+#include <vtkExtractCells.h>
+#include <vtkGeometryFilter.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkExtractEdges.h>
+#include <vtkImplicitDataSet.h>
+
+// debugging stuff
+#include <vtkDataSetMapper.h>
+#include <vtkSelectionNode.h>
+#include <vtkExtractSelection.h>
+#include <vtkSelection.h>
+#include <vtkVertexGlyphFilter.h>
+
+#define RENDERTEST
+//#define CONNECTTEST
+
+#ifdef RENDERTEST
+#include "vtkRectilinearGridGeometryFilter.h"
+#include "vtkRectilinearGridOutlineFilter.h"
+#include "vtkPolyDataMapper.h"
+#include "vtkActor.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRenderer.h"
+#include "vtkRenderWindow.h"
+#include "vtkCamera.h"
+#include "vtkProperty.h"
+#include <vtkDataSetMapper.h>
+#include <vtkSelectionNode.h>
+#include <vtkSelection.h>
+#include <vtkExtractSelection.h>
+#include <vtkExtractEdges.h>
+#include <vtkVertexGlyphFilter.h>
+#include <vtkTriangleFilter.h>
+#include <vtkImplicitHalo.h>
+#endif
+
+#ifdef CONNECTTEST
+#include <vtkVersion.h>
+#include <vtkSmartPointer.h>
+#include <vtkPolyData.h>
+#include <vtkCellData.h>
+#include <vtkDoubleArray.h>
+#include <vtkDataSet.h>
+#include <vtkSphereSource.h>
+#include <vtkTriangleFilter.h>
+#include <vtkExtractEdges.h>
+#include <vtkDataSetMapper.h>
+#include <vtkActor.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkSelectionNode.h>
+#include <vtkSelection.h>
+#include <vtkExtractSelection.h>
+#include <vtkProperty.h>
+#include <vtkVertexGlyphFilter.h>
+#endif
+
+vtkSmartPointer<vtkIdList> GetConnectedVertices(vtkSmartPointer<vtkDataSet> mesh, int id);
+vtkSmartPointer<vtkIdList> GetConnectedVertices2(vtkSmartPointer<vtkPolyData> mesh, int id);
+
+// This function is copyright (C) 2012 Joeseph Cappriotti
+vtkSmartPointer<vtkIdList> GetConnectedPoints(int id0, vtkSmartPointer<vtkDataSet> grid);
+
+using namespace Lemma;
+
+int main() {
+
+    // Define Sigma term
+    // Defaults to 1, (Poisson equation)
+
+    // Define source (G) term
+    //vtkSphere *Sphere = vtkSphere::New();
+    vtkImplicitHalo *Sphere = vtkImplicitHalo::New();
+        Sphere->SetCenter (0, 0, 0);
+        Sphere->SetRadius (60);
+        Sphere->SetFadeOut(.95);
+
+    ////////////////////////////////////////////
+    // Define solution mesh
+
+    // Create a rectilinear grid by defining three arrays specifying the
+    // coordinates in the x-y-z directions.
+
+    int    nx = 160;
+    double dx = 2.5;
+    double ox = -200;
+    vtkFloatArray *xCoords = vtkFloatArray::New();
+    for (int i=0; i<nx+1; i++) xCoords->InsertNextValue(ox + i*dx);
+
+    int    ny = 160;
+    double dy = 2.5;
+    double oy = -200;
+    vtkFloatArray *yCoords = vtkFloatArray::New();
+    for (int i=0; i<ny+1; i++) yCoords->InsertNextValue(oy + i*dy);
+
+    int    nz = 160;
+    double dz = 2.5;
+    double oz = -200;
+    vtkFloatArray *zCoords = vtkFloatArray::New();
+    for (int i=0; i<nz+1; i++) zCoords->InsertNextValue(oz + i*dz);
+
+    vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
+        rGrid->SetDimensions(nx+1, ny+1, nz+1);
+        rGrid->SetExtent(0, nx, 0, ny, 0,  nz);
+        rGrid->SetXCoordinates(xCoords);
+        rGrid->SetYCoordinates(yCoords);
+        rGrid->SetZCoordinates(zCoords);
+
+    vtkRectilinearGrid *sigGrid = vtkRectilinearGrid::New();
+        sigGrid->SetDimensions(nx+1, ny+1, nz+1);
+        sigGrid->SetExtent(0, nx, 0, ny, 0, nz);
+        sigGrid->SetXCoordinates(xCoords);
+        sigGrid->SetYCoordinates(yCoords);
+        sigGrid->SetZCoordinates(zCoords);
+        vtkDoubleArray *sigArray = vtkDoubleArray::New();
+        sigArray->SetNumberOfComponents(1);
+        for (int ii=0; ii<(nx+1)*(ny+1)*(nz+1); ++ii) {
+            sigArray->InsertTuple1(ii, 1.);
+        }
+        sigArray->SetName("sigma");
+        sigGrid->GetPointData()->AddArray(sigArray);
+        sigGrid->GetPointData()->SetScalars(sigArray);
+
+    vtkImplicitDataSet* implSigma = vtkImplicitDataSet::New();
+        implSigma->SetDataSet(sigGrid);
+        implSigma->SetOutValue(0.);
+        implSigma->SetOutGradient(0., 0., 0.);
+
+    FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
+        Solver->SetGFunction(Sphere);
+        //Solver->SetSigmaFunction(implSigma);
+        Solver->SetBoundaryStep(.5);
+        Solver->SetGrid(rGrid);
+        Solver->Solve("FEM_results.vtr");
+
+    vtkXMLDataSetWriter *Writer = vtkXMLDataSetWriter::New();
+        Writer->SetInputData(sigGrid);
+        std::string fname("sigma.vtr");
+        //Writer->Update();
+        //fname.append( Writer->GetDefaultFileExtension() );
+        Writer->SetFileName(fname.c_str());
+        Writer->Write();
+        Writer->Delete();
+
+    // Clean up
+    Sphere->Delete();
+    Solver->Delete();
+    rGrid->Delete();
+    xCoords->Delete();
+    yCoords->Delete();
+    zCoords->Delete();
+
+}
+
+// This function is copyright (C) 2012 Joeseph Cappriotti
+vtkSmartPointer<vtkIdList> GetConnectedPoints(int id0, vtkSmartPointer<vtkDataSet> grid){
+    vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
+    vtkSmartPointer<vtkIdList> cellList = vtkSmartPointer<vtkIdList>::New();
+    grid->GetPointCells(id0, cellList);
+    for(int i=0;i<cellList->GetNumberOfIds();++i){
+        vtkCell* cell = grid->GetCell(cellList->GetId(i));
+        if(cell->GetNumberOfEdges() > 0){
+            for(int j=0; j<cell->GetNumberOfEdges(); ++j){
+                vtkCell* edge = cell->GetEdge(j);
+                vtkIdList* edgePoints=edge->GetPointIds();
+                if(edgePoints->GetId(0)==id0){
+                    pointIds->InsertUniqueId(edgePoints->GetId(1));
+                }else if(edgePoints->GetId(1)==id0){
+                    pointIds->InsertUniqueId(edgePoints->GetId(0));
+                }
+            }
+        }
+    }
+    return pointIds;
+}
+
+vtkSmartPointer<vtkIdList> GetConnectedVertices(vtkSmartPointer<vtkDataSet> mesh, int id) {
+
+    std::cout << "number of points " << mesh->GetNumberOfPoints() << std::endl;
+    std::cout << "number of cells " << mesh->GetNumberOfCells() << std::endl;
+
+    vtkSmartPointer<vtkIdList> connectedVertices = vtkSmartPointer<vtkIdList>::New();
+
+    //get all cells that vertex 'id' is a part of
+    vtkSmartPointer<vtkIdList> cellIdList = vtkSmartPointer<vtkIdList>::New();
+    mesh->GetPointCells(id, cellIdList);   // cells using point
+    //mesh->GetCellPoints(id, cellIdList); // points defining cell
+
+    cout << "Vertex " << id << " is used in cell(s) ";
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
+        cout << cellIdList->GetId(i) << ", ";
+    }
+    cout << endl;
+
+    // TODO this is where things don't work
+    // for each cell, figure out what vertices are connected
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
+
+        cout << "\tcell id " << i << " : " << cellIdList->GetId(i) << endl;
+        std::cout << "\tcell has " << mesh->GetCell(i)->GetNumberOfPoints() << " points\n";
+
+        vtkExtractCells *Cell = vtkExtractCells::New();
+            Cell->SetInputData(mesh);
+            Cell->AddCellRange(i, i);
+            Cell->Update();
+
+        //std::cout << *Cell->GetOutput() << std::endl;
+
+        // extract surface
+        vtkDataSetSurfaceFilter *surf = vtkDataSetSurfaceFilter::New();
+            surf->UseStripsOn();
+            surf->SetInputData(Cell->GetOutput());
+            surf->Update();
+
+        //std::cout << *surf->GetOutput() << std::endl;
+
+        vtkSmartPointer<vtkTriangleFilter> triangleFilter =
+             vtkSmartPointer<vtkTriangleFilter>::New();
+        triangleFilter->SetInputData( surf->GetOutput() );
+            triangleFilter->Update();
+
+        vtkSmartPointer<vtkExtractEdges> extractEdges =
+            vtkSmartPointer<vtkExtractEdges>::New();
+        extractEdges->SetInputConnection(triangleFilter->GetOutputPort());
+            extractEdges->Update();
+
+        vtkSmartPointer<vtkPolyData> pmesh = extractEdges->GetOutput();
+        //vtkSmartPointer<vtkPolyData> pmesh = triangleFilter->GetOutput();
+        //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(surf->GetOutput(), id);
+//
+//         //vtkSmartPointer<vtkIdTypeArray> ids =
+//         //vtkSmartPointer<vtkIdTypeArray>::New();
+//         //ids->SetNumberOfComponents(1);
+//
+        vtkSmartPointer<vtkIdList> pointIdList = vtkSmartPointer<vtkIdList>::New();
+        //mesh->GetPointCells(cellIdList->GetId(i), pointIdList); // returns cells using point
+        mesh->GetCellPoints(cellIdList->GetId(i), pointIdList);   // points defining cell
+
+        //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(pmesh, id);
+        vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(surf->GetOutput(), id);
+        //vtkSmartPointer<vtkIdList> pconnectedVertices = GetConnectedVertices2(triangleFilter->GetOutput(), 1);
+        std::cout << "\tPoly Connected vertices: ";
+        for(vtkIdType ip = 0; ip < pconnectedVertices->GetNumberOfIds(); ip++) {
+            std::cout << pconnectedVertices->GetId(ip) << " ";
+            //ids->InsertNextValue(connectedVertices->GetId(i));
+            connectedVertices->InsertNextId( pointIdList->GetId( pconnectedVertices->GetId(ip) ) );
+        }
+        std::cout << std::endl;
+
+
+        cout << "\tcell " << i << " contains these points ";
+        for(vtkIdType ii = 0; ii < pointIdList->GetNumberOfIds(); ii++) {
+            cout << pointIdList->GetId(ii) << ", ";
+        }
+        cout << endl;
+
+        /*
+        cout << "\tEnd points are " << pointIdList->GetId(0) << " and " << pointIdList->GetId(1) << endl;
+
+        if(pointIdList->GetId(0) != id) {
+            //cout << "Connected to " << pointIdList->GetId(0) << endl;
+            connectedVertices->InsertNextId(pointIdList->GetId(0));
+        } else {
+            //cout << "Connected to " << pointIdList->GetId(1) << endl;
+            connectedVertices->InsertNextId(pointIdList->GetId(1));
+        }
+        */
+
+
+        #ifdef RENDERTEST
+        vtkSmartPointer<vtkDataSetMapper> sphereMapper =
+        vtkSmartPointer<vtkDataSetMapper>::New();
+            sphereMapper->SetInputConnection(extractEdges->GetOutputPort());
+        vtkSmartPointer<vtkActor> sphereActor =
+            vtkSmartPointer<vtkActor>::New();
+            sphereActor->SetMapper(sphereMapper);
+
+
+        vtkSmartPointer<vtkIdTypeArray> ids =
+            vtkSmartPointer<vtkIdTypeArray>::New();
+        ids->SetNumberOfComponents(1);
+
+        std::cout << "RENDERTEST Connected vertices: ";
+        for(vtkIdType i = 0; i < connectedVertices->GetNumberOfIds(); i++) {
+            std::cout << connectedVertices->GetId(i) << " ";
+            ids->InsertNextValue(connectedVertices->GetId(i));
+        }
+        std::cout << std::endl;
+
+        vtkSmartPointer<vtkDataSetMapper> connectedVertexMapper =
+        vtkSmartPointer<vtkDataSetMapper>::New();
+
+        {
+            vtkSmartPointer<vtkSelectionNode> selectionNode =
+                vtkSmartPointer<vtkSelectionNode>::New();
+            selectionNode->SetFieldType(vtkSelectionNode::POINT);
+            selectionNode->SetContentType(vtkSelectionNode::INDICES);
+            selectionNode->SetSelectionList(ids);
+
+            vtkSmartPointer<vtkSelection> selection =
+                vtkSmartPointer<vtkSelection>::New();
+            selection->AddNode(selectionNode);
+
+            vtkSmartPointer<vtkExtractSelection> extractSelection =
+                vtkSmartPointer<vtkExtractSelection>::New();
+
+            extractSelection->SetInputConnection(0, extractEdges->GetOutputPort());
+            extractSelection->SetInputData(1, selection);
+            extractSelection->Update();
+
+            vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
+                vtkSmartPointer<vtkVertexGlyphFilter>::New();
+            glyphFilter->SetInputConnection(extractSelection->GetOutputPort());
+            glyphFilter->Update();
+
+            connectedVertexMapper->SetInputConnection(glyphFilter->GetOutputPort());
+        }
+
+        vtkSmartPointer<vtkActor> connectedVertexActor =
+            vtkSmartPointer<vtkActor>::New();
+        connectedVertexActor->SetMapper(connectedVertexMapper);
+        connectedVertexActor->GetProperty()->SetColor(1,0,0);
+        connectedVertexActor->GetProperty()->SetPointSize(5);
+
+        vtkSmartPointer<vtkDataSetMapper> queryVertexMapper =
+            vtkSmartPointer<vtkDataSetMapper>::New();
+
+        {
+            vtkSmartPointer<vtkIdTypeArray> ids =
+                vtkSmartPointer<vtkIdTypeArray>::New();
+            ids->SetNumberOfComponents(1);
+            ids->InsertNextValue(id);
+
+            vtkSmartPointer<vtkSelectionNode> selectionNode =
+                vtkSmartPointer<vtkSelectionNode>::New();
+            selectionNode->SetFieldType(vtkSelectionNode::POINT);
+            selectionNode->SetContentType(vtkSelectionNode::INDICES);
+            selectionNode->SetSelectionList(ids);
+
+            vtkSmartPointer<vtkSelection> selection =
+                vtkSmartPointer<vtkSelection>::New();
+            selection->AddNode(selectionNode);
+
+            vtkSmartPointer<vtkExtractSelection> extractSelection =
+                vtkSmartPointer<vtkExtractSelection>::New();
+
+            extractSelection->SetInputConnection(0, extractEdges->GetOutputPort());
+            extractSelection->SetInputData(1, selection);
+            extractSelection->Update();
+
+            vtkSmartPointer<vtkVertexGlyphFilter> glyphFilter =
+                vtkSmartPointer<vtkVertexGlyphFilter>::New();
+            glyphFilter->SetInputConnection(extractSelection->GetOutputPort());
+            glyphFilter->Update();
+
+            queryVertexMapper->SetInputConnection(glyphFilter->GetOutputPort());
+        }
+
+        vtkSmartPointer<vtkActor> queryVertexActor =
+            vtkSmartPointer<vtkActor>::New();
+        queryVertexActor->SetMapper(queryVertexMapper);
+        queryVertexActor->GetProperty()->SetColor(0,1,0);
+        queryVertexActor->GetProperty()->SetPointSize(5);
+
+        // Create the usual rendering stuff.
+        vtkRenderer *renderer = vtkRenderer::New();
+        vtkRenderWindow *renWin = vtkRenderWindow::New();
+        renWin->AddRenderer(renderer);
+        vtkRenderWindowInteractor *iren = vtkRenderWindowInteractor::New();
+        iren->SetRenderWindow(renWin);
+
+        renderer->AddActor(queryVertexActor);
+        renderer->AddActor(connectedVertexActor);
+        renderer->AddActor(sphereActor);
+
+        //renderer->AddActor(wireActor);
+        //renderer->AddActor(owireActor);
+        renderer->SetBackground(.3,.2,.1);
+        renderer->ResetCamera();
+        renderer->GetActiveCamera()->Elevation(-150.0);
+        renderer->GetActiveCamera()->Azimuth(0.0);
+        renderer->GetActiveCamera()->Roll(0.0);
+        //renderer->GetActiveCamera()->Pitch(-45.0);
+        renderer->GetActiveCamera()->Zoom(1.0);
+
+        renWin->SetSize(300,300);
+
+        // interact with data
+        renWin->Render();
+        iren->Start();
+
+        // clean up rendering stuff
+        //plane->Delete();
+        //rgridMapper->Delete();
+        //wireActor->Delete();
+        //owireActor->Delete();
+        //outline->Delete();
+        //outlineMapper->Delete();
+        renderer->Delete();
+        renWin->Delete();
+        iren->Delete();
+        #endif // RENDERTEST
+
+        surf->Delete();
+        Cell->Delete();
+    }
+
+    return connectedVertices;
+}
+
+// from http://www.vtk.org/Wiki/VTK/Examples/Cxx/PolyData/VertexConnectivity
+vtkSmartPointer<vtkIdList> GetConnectedVertices2(vtkSmartPointer<vtkPolyData> mesh, int id) {
+
+    std::cout << "mesh points " << mesh->GetNumberOfPoints() << std::endl;
+
+    vtkSmartPointer<vtkIdList> connectedVertices =
+        vtkSmartPointer<vtkIdList>::New();
+
+    //get all cells that vertex 'id' is a part of
+    vtkSmartPointer<vtkIdList> cellIdList =
+        vtkSmartPointer<vtkIdList>::New();
+    mesh->GetPointCells(id, cellIdList);
+
+    /*
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
+        cout << cellIdList->GetId(i) << ", ";
+    }
+    cout << endl;
+    */
+
+    for(vtkIdType i = 0; i < cellIdList->GetNumberOfIds(); i++) {
+        //cout << "id " << i << " : " << cellIdList->GetId(i) << endl;
+
+        vtkSmartPointer<vtkIdList> pointIdList =
+            vtkSmartPointer<vtkIdList>::New();
+        mesh->GetCellPoints(cellIdList->GetId(i), pointIdList);
+        //cout << "End points are " << pointIdList->GetId(0) << " and " << pointIdList->GetId(1) << endl;
+        if(pointIdList->GetId(0) != id) {
+            //cout << "Connected to " << pointIdList->GetId(0) << endl;
+            connectedVertices->InsertNextId(pointIdList->GetId(0));
+        } else {
+            //cout << "Connected to " << pointIdList->GetId(1) << endl;
+            connectedVertices->InsertNextId(pointIdList->GetId(1));
+        }
+    }
+    return connectedVertices;
+}

examples/utFEM4EllipticPDE_bhmag.cpp

@@ -0,0 +1,173 @@
+// ===========================================================================
+//
+//       Filename:  utFEM4EllipticPDE.cpp
+//
+//        Created:  08/16/12 19:49:10
+//       Compiler:  Tested with g++, icpc, and MSVC 2010
+//
+//         Author:  Trevor Irons (ti)
+//
+//   Organisation:  Colorado School of Mines (CSM)
+//                  United States Geological Survey (USGS)
+//
+//          Email:  tirons@mines.edu, tirons@usgs.gov
+//
+//  This program is free software: you can redistribute it and/or modify
+//  it under the terms of the GNU General Public License as published by
+//  the Free Software Foundation, either version 3 of the License, or
+//  (at your option) any later version.
+//
+//  This program is distributed in the hope that it will be useful,
+//  but WITHOUT ANY WARRANTY; without even the implied warranty of
+//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+//  GNU General Public License for more details.
+//
+//  You should have received a copy of the GNU General Public License
+//  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// ===========================================================================
+
+/**
+  @file
+  @author   Trevor Irons
+  @date     08/16/12
+  @version   0.0
+ **/
+
+#include "FEM4EllipticPDE.h"
+#include "vtkSphere.h"
+#include "vtkRectilinearGrid.h"
+#include "vtkFloatArray.h"
+#include <vtkIdList.h>
+#include <vtkIdTypeArray.h>
+#include <vtkCell.h>
+#include <vtkTriangleFilter.h>
+#include <vtkDataSetSurfaceFilter.h>
+#include <vtkExtractCells.h>
+#include <vtkGeometryFilter.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkExtractEdges.h>
+#include <vtkImplicitDataSet.h>
+#include <vtkXMLUnstructuredGridReader.h>
+
+// debugging stuff
+#include <vtkDataSetMapper.h>
+#include <vtkSelectionNode.h>
+#include <vtkExtractSelection.h>
+#include <vtkSelection.h>
+#include <vtkVertexGlyphFilter.h>
+
+#define RENDERTEST
+//#define CONNECTTEST
+
+#ifdef RENDERTEST
+#include "vtkRectilinearGridGeometryFilter.h"
+#include "vtkRectilinearGridOutlineFilter.h"
+#include "vtkPolyDataMapper.h"
+#include "vtkActor.h"
+#include "vtkRenderWindowInteractor.h"
+#include "vtkRenderer.h"
+#include "vtkRenderWindow.h"
+#include "vtkCamera.h"
+#include "vtkProperty.h"
+#include <vtkDataSetMapper.h>
+#include <vtkSelectionNode.h>
+#include <vtkSelection.h>
+#include <vtkExtractSelection.h>
+#include <vtkExtractEdges.h>
+#include <vtkVertexGlyphFilter.h>
+#include <vtkTriangleFilter.h>
+#include <vtkImplicitHalo.h>
+#endif
+
+using namespace Lemma;
+
+Real Magnet(const Real& x, const Real& y, const Real& z);
+
+int main(int argc, char**argv) {
+
+    if (argc < 4) {
+        std::cout << "usage: ./utFEMEllipticPDE_bhmag  <g.vtu>  <mesh.vtu>  <results.vtu>" << std::endl;
+        //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
+    double ox = -.175;
+    vtkFloatArray *xCoords = vtkFloatArray::New();
+    for (int i=0; i<nx+1; i++) xCoords->InsertNextValue(ox + i*dx);
+
+    int    ny = 80;
+    double dy =  .004375;
+    double oy = -.175;
+    vtkFloatArray *yCoords = vtkFloatArray::New();
+    for (int i=0; i<ny+1; i++) yCoords->InsertNextValue(oy + i*dy);
+
+    int    nz = 343; // 685; // 160
+    double dz = .004375;
+    double oz = 9.75;
+    vtkFloatArray *zCoords = vtkFloatArray::New();
+    for (int i=0; i<nz+1; i++) zCoords->InsertNextValue(oz + i*dz);
+
+    vtkRectilinearGrid *rGrid = vtkRectilinearGrid::New();
+        rGrid->SetDimensions(nx+1, ny+1, nz+1);
+        rGrid->SetExtent(0, nx, 0, ny, 0,  nz);
+        rGrid->SetXCoordinates(xCoords);
+        rGrid->SetYCoordinates(yCoords);
+        rGrid->SetZCoordinates(zCoords);
+
+    ////////////////////////////////////////////
+    // Define Sigma/mu term
+    // TODO
+
+    ////////////////////////////////////////////
+    // Define source (G) term
+    vtkXMLUnstructuredGridReader* GReader = vtkXMLUnstructuredGridReader::New();
+        std::cout << "Reading G file " << argv[1] << std::endl;
+        GReader->SetFileName(argv[1]);
+        GReader->Update();
+
+    vtkImplicitDataSet* implG = vtkImplicitDataSet::New();
+        implG->SetDataSet(GReader->GetOutput());
+        implG->SetOutValue(0.);
+        implG->SetOutGradient(0., 0., 0.);
+
+    ////////////////////////////////////////////
+    // Define solution mesh
+    vtkXMLUnstructuredGridReader* MeshReader = vtkXMLUnstructuredGridReader::New();
+        std::cout << "Reading Mesh file " << argv[2] << std::endl;
+        MeshReader->SetFileName(argv[2]);
+        MeshReader->Update();
+
+    ////////////////////////////////////////////
+    // Solve
+    FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
+        Solver->SetGFunction(implG);
+        //Solver->SetGFunction(Magnet);
+        //Solver->SetSigmaFunction(implSigma);
+        Solver->SetBoundaryStep(.05);
+        Solver->SetGrid(MeshReader->GetOutput());
+        //Solver->SetGrid(rGrid);
+        Solver->Solve(argv[3]);
+
+    // Clean up
+    Solver->Delete();
+    GReader->Delete();
+    implG->Delete();
+
+    exit(EXIT_SUCCESS);
+}
+
+Real Magnet(const Real& x, const Real& y, const Real& z) {
+
+    if (z < 10 || z > 11 ) {
+        return 0.;
+    } else if ( std::sqrt(x*x + y*y) <= 0.05 ) {
+        if (x>0.) return  1.;
+        if (x<0.) return -1.;
+    }
+    return 0.;
+}
+

examples/utVTKDC.cpp

@@ -0,0 +1,104 @@
+/* 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 <cmath>
+#include <Eigen/Core>
+
+#include "DCSurvey.h"
+#include "FEM4EllipticPDE.h"
+
+#ifdef HAVE_YAMLCPP
+#include "yaml-cpp/yaml.h"
+#endif
+
+using namespace Lemma;
+
+int main(int argc, char**argv) {
+
+    std::cout << "Mesh processing routine\n";
+    if (argc<4) {
+        std::cout << "usage:\n" << "./utVTKEdge  mesh.vtu  input.yaml   results.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";
+
+    std::ifstream ifstr(argv[2]);
+    //std::vector<YAML::Node> nodes = YAML::LoadAll(ifstr);
+    YAML::Node nodes = YAML::Load(ifstr);
+    DCSurvey* Survey = DCSurvey::DeSerialize(nodes);
+    //std::cout << *Survey << std::endl;
+
+
+    ////////////////////////////////////////////
+    // Solve
+    FEM4EllipticPDE *Solver = FEM4EllipticPDE::New();
+        //Solver->SetGFunction(implG);
+        //Solver->SetGFunction(Magnet);
+        //Solver->SetSigmaFunction(implSigma);
+        //Solver->SetBoundaryStep(.05);
+
+        Solver->SetGrid(uGrid);
+        Solver->SetupDC(Survey, 0);
+        Solver->Solve( argv[3] );
+        //Solver->SetGrid(rGrid);
+        //Solver->Solve(argv[3]);
+
+    Survey->Delete();
+    Solver->Delete();
+
+    exit(EXIT_SUCCESS);
+}

examples/utVTKEdgeG.cpp

@@ -0,0 +1,95 @@
+/* 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 <vtkXMLUnstructuredGridReader.h>
+#include <vtkXMLUnstructuredGridWriter.h>
+#include <vtkDoubleArray.h>
+#include <vtkPointData.h>
+
+#include <cmath>
+
+
+int main(int argc, char**argv) {
+
+    std::cout << "Mesh processing routine\n";
+    if (argc<3) {
+        std::cout << "usage:\n" << "./utVTKEdge  inMesh.vtu outG.vtu" << std::endl;
+    }
+
+    vtkXMLUnstructuredGridReader* Reader = vtkXMLUnstructuredGridReader::New();
+        std::cout << "Reading" << argv[1] << std::endl;
+        Reader->SetFileName(argv[1]);
+        Reader->Update();
+
+        int nc = Reader->GetOutput()->GetNumberOfCells()  ;
+        int nn = Reader->GetOutput()->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 = .25;          // radius of magnet
+    double D0 = 10;          // Top of magnet
+    double D1 = 11;          // Bottom of magnet
+    double eps = 1e-2;       // epsilon for on the point
+    // Pol = double[3];
+
+    // Declare new array for the data
+    vtkDoubleArray* G = vtkDoubleArray::New();
+        G->SetNumberOfComponents(1);
+        G->SetNumberOfTuples(nn);
+        G->SetName("G");
+
+    // Loop over nodes, look at position, set with G if on boundary
+    double point[3];
+    for (int in=0; in<nn; ++in) {
+        Reader->GetOutput()->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] );
+        if ( raddist > R - eps && raddist < R + eps && point[2] < D1 && point[2] > D0 ) {
+            if ( point[1] < 0 ) {
+                G->InsertTuple1( in,   1 );
+            } else if (point[1] > 0) {
+                G->InsertTuple1( in,  -1 );
+            } else {
+                G->InsertTuple1( in, 0 );
+            }
+        } else {
+            G->InsertTuple1( in, 0 );
+        }
+
+    }
+
+    // Add new data
+    Reader->GetOutput()->GetPointData()->SetScalars( G );
+
+    // Write out new file
+    vtkXMLUnstructuredGridWriter* Writer = vtkXMLUnstructuredGridWriter::New();
+        Writer->SetInputConnection(Reader->GetOutputPort());
+        Writer->SetFileName( argv[2] );
+        Writer->Write();
+
+    Reader->Delete();
+
+
+}

examples/utVTKEdgeGsphere.cpp

@@ -0,0 +1,172 @@
+/* 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 <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(nn);
+        G->SetName("G");
+
+    vtkDoubleArray* phi = vtkDoubleArray::New();
+        phi->SetNumberOfComponents(1);
+        phi->SetNumberOfTuples(nn);
+        phi->SetName("analytic_phi");
+
+    // Loop over nodes, look at position, set with G if on boundary
+    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
+        if ( std::string(argv[3]) == "gravity") {
+            if ( raddist < R + eps) {
+                G->InsertTuple1( in, 1 );
+            } else {
+                G->InsertTuple1( in, 0 );
+            }
+        }
+
+        else 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 );
+            }
+        }
+
+        // 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);
+            if (raddist < R) {
+                phi->InsertTuple1( in, mass * ( 3*(R*R) - raddist*raddist )/(2*(R*R*R)) ); // (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 );
+            }
+        }
+    }
+
+    // Add new data
+    uGrid->GetPointData()->AddArray( phi );
+    uGrid->GetPointData()->SetScalars( 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/utmerge.cpp

@@ -0,0 +1,126 @@
+/* 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      07/29/2014 02:01:21 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 <vtkVersion.h>
+#include <vtkSmartPointer.h>
+#include <vtkPolyData.h>
+#include <vtkSphereSource.h>
+#include <vtkConeSource.h>
+#include <vtkXMLPolyDataReader.h>
+#include <vtkCleanPolyData.h>
+#include <vtkAppendPolyData.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkActor.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderer.h>
+#include <vtkRenderWindowInteractor.h>
+
+int main(int argc, char *argv[])
+{
+  vtkSmartPointer<vtkPolyData> input1 =
+    vtkSmartPointer<vtkPolyData>::New();
+  vtkSmartPointer<vtkPolyData> input2 =
+    vtkSmartPointer<vtkPolyData>::New();
+
+  if(argc == 1) //command line arguments not specified
+    {
+    vtkSmartPointer<vtkSphereSource> sphereSource =
+      vtkSmartPointer<vtkSphereSource>::New();
+    sphereSource->SetCenter(5,0,0);
+    sphereSource->Update();
+
+    input1->ShallowCopy(sphereSource->GetOutput());
+
+    vtkSmartPointer<vtkConeSource> coneSource =
+      vtkSmartPointer<vtkConeSource>::New();
+    coneSource->Update();
+
+    input2->ShallowCopy(coneSource->GetOutput());
+    }
+  else
+  {
+    if(argc != 3)
+      {
+      std::cout << "argc = " << argc << std::endl;
+      std::cout << "Required arguments: File1 File2" << std::endl;
+      return EXIT_FAILURE;
+      }
+    std::string inputFilename1 = argv[1];
+    std::string inputFilename2 = argv[2];
+    vtkSmartPointer<vtkXMLPolyDataReader> reader1 =
+      vtkSmartPointer<vtkXMLPolyDataReader>::New();
+    reader1->SetFileName(inputFilename1.c_str());
+    reader1->Update();
+    input1->ShallowCopy(reader1->GetOutput());
+
+    vtkSmartPointer<vtkXMLPolyDataReader> reader2 =
+      vtkSmartPointer<vtkXMLPolyDataReader>::New();
+    reader2->SetFileName(inputFilename2.c_str());
+    reader2->Update();
+    input2->ShallowCopy(reader2->GetOutput());
+  }
+
+  //Append the two meshes
+  vtkSmartPointer<vtkAppendPolyData> appendFilter =
+    vtkSmartPointer<vtkAppendPolyData>::New();
+#if VTK_MAJOR_VERSION <= 5
+  appendFilter->AddInputConnection(input1->GetProducerPort());
+  appendFilter->AddInputConnection(input2->GetProducerPort());
+#else
+  appendFilter->AddInputData(input1);
+  appendFilter->AddInputData(input2);
+#endif
+  appendFilter->Update();
+
+  // Remove any duplicate points.
+  vtkSmartPointer<vtkCleanPolyData> cleanFilter =
+    vtkSmartPointer<vtkCleanPolyData>::New();
+  cleanFilter->SetInputConnection(appendFilter->GetOutputPort());
+  cleanFilter->Update();
+
+  //Create a mapper and actor
+  vtkSmartPointer<vtkPolyDataMapper> mapper =
+    vtkSmartPointer<vtkPolyDataMapper>::New();
+  mapper->SetInputConnection(cleanFilter->GetOutputPort());
+
+  vtkSmartPointer<vtkActor> actor =
+    vtkSmartPointer<vtkActor>::New();
+  actor->SetMapper(mapper);
+
+  //Create a renderer, render window, and interactor
+  vtkSmartPointer<vtkRenderer> renderer =
+    vtkSmartPointer<vtkRenderer>::New();
+  vtkSmartPointer<vtkRenderWindow> renderWindow =
+    vtkSmartPointer<vtkRenderWindow>::New();
+  renderWindow->AddRenderer(renderer);
+  vtkSmartPointer<vtkRenderWindowInteractor> renderWindowInteractor =
+    vtkSmartPointer<vtkRenderWindowInteractor>::New();
+  renderWindowInteractor->SetRenderWindow(renderWindow);
+
+  //Add the actors to the scene
+  renderer->AddActor(actor);
+  renderer->SetBackground(.3, .2, .1); // Background color dark red
+
+  //Render and interact
+  renderWindow->Render();
+  renderWindowInteractor->Start();
+
+  return EXIT_SUCCESS;
+}
+