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 <vtkCellData.h>
#include <vtkCell.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");
        if ( std::string(argv[3]) == "mag") {
            G->SetName("kappa");
            G->SetNumberOfTuples(nc);
        }
        else {
            G->SetName("G");
            G->SetNumberOfTuples(nn);
        }
    vtkDoubleArray* phi = vtkDoubleArray::New();
        phi->SetNumberOfComponents(1);
        phi->SetNumberOfTuples(nn);
        phi->SetName("analytic_phi");

    double point[3];


    if ( std::string(argv[3]) == "mag") {
        // Loop over cells and decide if all of them lie within sphere
        for (int ic=0; ic<nc; ++ic) {
            assert ( uGrid->GetCell(ic)->GetNumberOfPoints() == 4 );
            // TODO, in production code we might not want to do this check here
            if ( uGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
                throw std::runtime_error("Non-tetrahedral mesh encountered!");
            }
            // construct coordinate matrix C
            //Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
            bool cellin(true);
            for (int ii=0; ii<4; ++ii) {
                double* point =  uGrid->GetCell(ic)->GetPoints()->GetPoint(ii); //[ipc] ;
                if ( sqrt( point[0]*point[0] + point[1]*point[1] + point[2]*point[2] ) > R+eps ) {
                    cellin = false;
                }
            }
            if (cellin) {
                G->InsertTuple1( ic, 1 );
            } else {
                G->InsertTuple1( ic, 0 );
            }
        }
        uGrid->GetCellData()->AddArray( G );

    } else {
        // Loop over nodes, look at position, set with G if on boundary
        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 );
                }
            }
            if ( std::string(argv[3]) == "magnet") {
                if ( raddist > R - eps && raddist < R + eps ) {
                    // \rho = \nabla \cdot \mathbf{M}
                    // Use divergence theorm --> \mahtbf{M} \cdot \hat{n}
                    Eigen::Vector3d n;
                    n << point[0],point[1],point[2];
                    n.array() /= raddist;
                    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] );
                    //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 );
                }
            }
        }

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


}