src/LinearFieldInterpolation.cpp

Go to the documentation of this file.

//
// C++ Implementation: %{MODULE}
//
// Description:
//
//
// Author: %{AUTHOR} <%{EMAIL}>, (C) %{YEAR}
//
// Copyright: See COPYING file that comes with this distribution
//
//

using namespace std;

#include "myrlog.h"
#include "LinearFieldInterpolation.h"



namespace postprocess {

  //using namespace std;
using namespace mesh;

LinearFieldInterpolation::LinearFieldInterpolation(const NedelecMesh& nedelec_mesh, const colarray::Vector<double>& q)
    : nedelec_mesh_(nedelec_mesh)
{
    TetMesh* mesh(nedelec_mesh_.getTetMesh());
    vector<int> node_count_(mesh->get_nof_points(), 0);

    // Initialise field values with zero vectors    
    node_values_.resize(mesh->get_nof_points(), Vector3(0.0, 0.0, 0.0));
     
    // Evaluate field at mesh points    
    for (TetMesh::tet_iterator it = mesh->tet_begin(); it != mesh->tet_end(); ++ it) {
        Tet& tet = *it;
        NedelecElement* elem = nedelec_mesh_.get_element(tet.get_id());
        
        // build vector of local DOFs
        const int* dof_ids = elem->get_dof_ids(mesh, &tet);
        colarray::Vector<double> dofs(elem->get_nof_dofs());
        for (int i = 0; i < elem->get_nof_dofs(); i ++) {
            int mapped_id = nedelec_mesh_.map_dof(dof_ids[i]);
            if (mapped_id == -1)
                dofs(i) = 0.0;
            else
                dofs(i) = q(mapped_id);
        }

        for (int i = 0; i < 4; ++ i) {
            mesh::id_t point_id = tet.get_corner_id(i);
            node_values_[point_id] += elem->eval(&tet, mesh->get_point(point_id)->get_coord(), dofs);
            node_count_[point_id] += 1;
        }
    }
    
    // Build average
    for (mesh::id_t i = 0; i < mesh->get_nof_points(); ++ i) {
        node_values_[i] /= node_count_[i];
    }
}

LinearFieldInterpolation::~LinearFieldInterpolation()
{
}

mesh::Vector3 postprocess::LinearFieldInterpolation::eval(const mesh::Vector3& x)
{
    // identify tetrahedr(on|a) containing x using octree data structure
    TetMesh* tet_mesh(nedelec_mesh_.getTetMesh());
    std::vector<mesh::Tet *> tets;
    tet_mesh->find_tets_by_point(x, tets);
    // warn if no x is not in mesh
    if (tets.size() == 0) {
        rError("eval: No tet found at location (%g,%g,%g).", x.x, x.y, x.z);
        return Vector3::ZERO;
    }
    // Choose first tet. The choice is not critical here, since the interpolated 
    // field is contignous
    mesh::Tet *tet(tets[0]);
    // Compute simplex coordinates in chosen Tet.
    Vector4 x_simplex;
    tet->cartesian_to_simplex(x, x_simplex);
    // Sum field contributions from all 4 corners: field = sum(node_values_[point_id[i]] * x_simplex[i], i==0..3)
    Vector3 field(Vector3::ZERO);
    for (int i = 0; i < 4; ++ i)
        field += node_values_[tet->get_corner_id(i)] * x_simplex[i];
    return field;
}

mesh::Vector3 postprocess::LinearFieldInterpolation::eval_curl(const mesh::Vector3& x)
{
    rError("LinearFieldInterpolation::eval_curl not implemented.");
    return Vector3::ZERO;
}

} // namespace postprocess

---