FEM4EllipticPDE/examples/FEM4EllipticPDE_bhmag.cpp

// ===========================================================================
//
//       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) {

    Eigen::initParallel();

    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->SetupPotential();
        //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.;
}