FEM4EllipticPDE/src/LinearMag.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      03/21/2016 02:10:08 PM
 * @author    Trevor Irons (ti)
 * @email     tirons@egi.utah.edu
 * @copyright Copyright (c) 2016, University of Utah
 * @copyright Copyright (c) 2016, Lemma Software, LLC
 */

#include "LinearMag.h"

namespace Lemma {

// ====================  FRIEND METHODS  =====================
std::ostream &operator << (std::ostream &stream, const LinearMag &ob) {
    stream << ob.Serialize()  << "\n";
    return stream;
}

// ====================  LIFECYCLE     =======================

//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  LinearMag
// Description:  constructor (protected)
//--------------------------------------------------------------------------------------
LinearMag::LinearMag (const ctor_key& key) : FEM4EllipticPDE(key) {

}  // -----  end of method LinearMag::LinearMag  (constructor)  -----

//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  LinearMag
// Description:  DeSerializing constructor (protected)
//--------------------------------------------------------------------------------------
LinearMag::LinearMag (const YAML::Node& node, const ctor_key& key) : FEM4EllipticPDE(node, key) {

}  // -----  end of method LinearMag::LinearMag  (constructor)  -----

//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  New()
// Description:  public constructor
//--------------------------------------------------------------------------------------
std::shared_ptr<LinearMag> LinearMag::NewSP() {
    return std::make_shared<LinearMag>( ctor_key() );
}

//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  ~LinearMag
// Description:  destructor (protected)
//--------------------------------------------------------------------------------------
LinearMag::~LinearMag () {

}  // -----  end of method LinearMag::~LinearMag  (destructor)  -----

//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  Serialize
//--------------------------------------------------------------------------------------
YAML::Node  LinearMag::Serialize (  ) const {
    YAML::Node node = FEM4EllipticPDE::Serialize();;
    node.SetTag( this->GetName() );
    // FILL IN CLASS SPECIFICS HERE
    node["B0"] = B0;
    return node;
}		// -----  end of method LinearMag::Serialize  -----


//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  DeSerialize
//--------------------------------------------------------------------------------------
std::shared_ptr<LinearMag> LinearMag::DeSerialize ( const YAML::Node& node ) {
    if ( node.Tag() != "LinearMag") {
        throw DeSerializeTypeMismatch( "LinearMag", node.Tag() );
    }
    return std::make_shared<LinearMag>( node, ctor_key() );
}		// -----  end of method LinearMag::DeSerialize  -----

//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  SetInducingMagField
//--------------------------------------------------------------------------------------
void LinearMag::SetInducingMagField ( const Real& intensity, const Real& inc,
                const Real& dec, const MAGUNITS& U ) {

    B0(0) = intensity * std::cos(inc) * std::cos(dec); // northing
    B0(1) = intensity * std::cos(inc) * std::sin(dec); // easting
    B0(2) = intensity * std::sin(inc);                 // z
    ScaleB0(U);
    return ;
}		// -----  end of method LinearMag::SetInducingMagField  -----


//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  SetInducingMagFieldVector
//--------------------------------------------------------------------------------------
void LinearMag::SetInducingMagFieldVector ( const Vector3r& BB0, const MAGUNITS& U  ) {
    B0 = BB0;
    ScaleB0(U);
    return ;
}		// -----  end of method LinearMag::SetInducingMagFieldVector  -----


//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  ScaleB0
//--------------------------------------------------------------------------------------
void LinearMag::ScaleB0 ( const MAGUNITS& U ) {
    switch ( U ) {
        case TESLA:
            break;
        case NANOTESLA:
            B0 *= 1e-9;
            break;
        case GAUSS:
            B0 *= 1e-4;
            break;
    }
    return ;
}		// -----  end of method LinearMag::ScaleB0  -----


//--------------------------------------------------------------------------------------
//       Class:  LinearMag
//      Method:  CalculateRHS
//--------------------------------------------------------------------------------------
void LinearMag::CalculateRHS ( const std::string& susName ) {

    std::cout << "Calculating RHS...";
    std::cout.flush();

    if ( !vtkGrid->GetCellData()->GetScalars(susName.c_str()) ) {
        std::string err("No cell data by name ");
        err.append(susName);
        throw std::runtime_error(err.c_str());
    }

    if (!vtkGrid->GetNumberOfPoints()) {
        throw std::runtime_error("Number of points zero in input grid!");
    }

    vtkDoubleArray* G = vtkDoubleArray::New();
    G->SetNumberOfComponents(1);
    G->SetNumberOfTuples( vtkGrid->GetNumberOfPoints() );
    G->SetName("G");
    //g.resize(vtkGrid->GetNumberOfPoints());
    VectorXr GG = VectorXr::Zero( vtkGrid->GetNumberOfPoints() );

    // Iterate over all the points or all of the cells?
    for (int ic=0; ic < vtkGrid->GetNumberOfCells(); ++ic) {

        if ( vtkGrid->GetCell(ic)->GetNumberOfPoints() != 4 ) {
            throw std::runtime_error("Non-tetrahedral mesh encountered!");
        }
        Real cellSus = vtkGrid->GetCellData()->GetScalars(susName.c_str())->GetTuple(ic)[0];

        Eigen::Matrix<Real, 4, 4> C = Eigen::Matrix<Real, 4, 4>::Zero() ;
        for (int ii=0; ii<4; ++ii) {
            double* pts =  vtkGrid->GetCell(ic)->GetPoints()->GetPoint(ii);
            C(ii, 0) = 1;
            C(ii, 1) =  pts[0];
            C(ii, 2) =  pts[1];
            C(ii, 3) =  pts[2];
        }

        /* The indices */
        vtkIdList* Ids = vtkGrid->GetCell(ic)->GetPointIds();
        int ID[4];
            ID[0] = Ids->GetId(0);
            ID[1] = Ids->GetId(1);
            ID[2] = Ids->GetId(2);
            ID[3] = Ids->GetId(3);

        /* the 4 faces of the tetrahedra
            ID[0] ID[1] ID[2]
            ID[0] ID[1] ID[3]
            ID[0] ID[2] ID[3]
            ID[1] ID[2] ID[3]
        */
        // Face 0, ID 0,1,2
        Eigen::Matrix<Real, 3, 2> CC = Eigen::Matrix<Real, 3, 2>::Ones() ;
            CC.col(1) = C.row(0).tail<3>() - C.row(1).tail<3>();
            CC.col(2) = C.row(0).tail<3>() - C.row(2).tail<3>();
            Vector3r nhat = CC.col(1).cross(CC.col(2));
            nhat.array() /= nhat.norm();
            Real flux = cellSus*nhat.dot(B0);
            GG(ID[0]) += flux;
            GG(ID[1]) += flux;
            GG(ID[2]) += flux;
	// Face 1, ID 0,1,3
        {
            CC.col(1) = C.row(0).tail<3>() - C.row(1).tail<3>();
            CC.col(2) = C.row(0).tail<3>() - C.row(3).tail<3>();
            Vector3r nhat = CC.col(1).cross(CC.col(2));
            nhat.array() /= nhat.norm();
            Real flux = cellSus*nhat.dot(B0);
            GG(ID[0]) += flux;
            GG(ID[1]) += flux;
            GG(ID[3]) += flux;
        }
        // Face 2, ID 0,2,3
        {
            CC.col(1) = C.row(0).tail<3>() - C.row(2).tail<3>();
            CC.col(2) = C.row(0).tail<3>() - C.row(3).tail<3>();
            Vector3r nhat = CC.col(1).cross(CC.col(2));
            nhat.array() /= nhat.norm();
            Real flux = cellSus*nhat.dot(B0);
            GG(ID[0]) += flux;
            GG(ID[2]) += flux;
            GG(ID[3]) += flux;
        }
        // Face 3, ID 1,2,3
        {
            CC.col(1) = C.row(1).tail<3>() - C.row(2).tail<3>();
            CC.col(2) = C.row(1).tail<3>() - C.row(3).tail<3>();
            Vector3r nhat = CC.col(1).cross(CC.col(2));
            nhat.array() /= nhat.norm();
            Real flux = cellSus*nhat.dot(B0);
            GG(ID[1]) += flux;
            GG(ID[2]) += flux;
            GG(ID[3]) += flux;
        }

    }
    for (int ip=0; ip<vtkGrid->GetNumberOfPoints(); ++ip) {
        G->InsertTuple1( ip, GG(ip) );
    }

    vtkGrid->GetPointData()->AddArray( G );
    vtkGrid->GetPointData()->SetScalars( G );
    std::cout << "finished" << std::endl;

    return ;
}		// -----  end of method LinearMag::CalculateRHS  -----


}		// -----  end of Lemma  name  -----