FEM4EllipticPDE/examples/VTKEdgeGsphere.cpp

/* 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); // volume * density (1)
            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();


}