faceted_setup.h

Go to the documentation of this file.

/*
  This file is part of the Ristra portage project.
  Please see the license file at the root of this repository, or at:
  https://github.com/laristra/portage/blob/master/LICENSE
*/

#ifndef FACETED_SETUP_H
#define FACETED_SETUP_H

#include <vector>
#include <algorithm>
#include <stdexcept>

#include "wonton/mesh/AuxMeshTopology.h"
#include "wonton/support/Point.h"
#include "portage/support/weight.h"
#include "portage/support/portage.h"

namespace Portage { namespace Meshfree { namespace Weight {

template<int dim, class Mesh>
void faceted_setup_cell(Mesh const& mesh,
                        Portage::vector<std::vector<std::vector<double>>>& smoothing_lengths,
                        Portage::vector<Wonton::Point<dim>>& extents,
                        double smoothing_factor, double boundary_factor) {

  if (smoothing_factor < 0.) {
    throw std::runtime_error("smoothing_factor must be non-negative");
  } else if (boundary_factor < 0.) {
    throw std::runtime_error("boundary_factor must be non-negative");
  }

  int ncells = mesh.num_owned_cells();
  smoothing_lengths.resize(ncells);
  extents.resize(ncells);

  // get face normal and distance information
  for (int i = 0; i < ncells; i++) {
    std::vector<int> faces, fdirs;
    mesh.cell_get_faces_and_dirs(i, &faces, &fdirs);

    int nfaces = faces.size();
    std::vector<std::vector<double>> h = smoothing_lengths[i];
    h.resize(nfaces);

    Wonton::Point<dim> fcent, ccent, normal, distance;
    std::vector<int> fnodes;
    std::vector<Wonton::Point<dim>> fncoord;
    for (int j = 0; j < nfaces; j++) {
      // get face data
      mesh.face_centroid(faces[j], &fcent);
      mesh.cell_centroid(i, &ccent);
      mesh.face_get_nodes(faces[j], &fnodes);
      int const num_face_nodes = fnodes.size();
      fncoord.resize(num_face_nodes + 1);
      for (int k = 0; k < num_face_nodes; k++) {
        mesh.node_get_coordinates(fnodes[k], &(fncoord[k]));
      }
      fncoord[num_face_nodes] = fncoord[0];

      if (dim == 2) {
        // Get 2d face normal
        std::vector<double> delta(2);
        for (int k = 0; k < 2; k++) delta[k] = fncoord[1][k] - fncoord[0][k];
        normal[0] = delta[1];
        normal[1] = -delta[0];
        if (fdirs[j] < 0) normal *= -1.0;
        double norm = sqrt(normal[0] * normal[0] + normal[1] * normal[1]);
        normal /= norm;
        h[j].resize(3);
        h[j][0] = normal[0];
        h[j][1] = normal[1];
      } else if (dim == 3) {
        // Get 3d face normal using average cross product, in case face is not flat
        // or has nodes close together.
        std::vector<std::vector<double>> bivec(2, std::vector<double>(3, 0.));
        for (int k = 0; k < 3; k++) {
          bivec[0][k] = fncoord[0][k] - fcent[k];
          normal[k] = 0.;
        }

        //int const num_face_nodes = fnodes.size();
        for (int k = 1; k <= num_face_nodes; k++) {
          for (int m = 0; m < 3; m++) bivec[1][m] = fncoord[k][m] - fcent[m];
          std::vector<double> cross(3, 0.);
          cross[0] = bivec[1][2] * bivec[0][1] - bivec[1][1] * bivec[0][2];
          cross[1] = -bivec[1][2] * bivec[0][0] + bivec[1][0] * bivec[0][2];
          cross[2] = bivec[1][1] * bivec[0][0] - bivec[1][0] * bivec[0][1];
          double norm = sqrt(cross[0] * cross[0] + cross[1] * cross[1] + cross[2] * cross[2]);
          for (int m = 0; m < 3; m++) {
            cross[m] /= norm;
          }
          normal += cross;
          bivec[0] = bivec[1];
        }
        if (fdirs[j] < 0) normal *= -1.0;
        double norm = 0.;
        for (int k = 0; k < 3; k++) norm += normal[k] * normal[k];
        norm = sqrt(norm);
        normal /= norm;
        h[j].resize(4);
        for (int k = 0; k < 3; k++) h[j][k] = normal[k];
      }

      distance = Wonton::Point<dim>(fcent - ccent);
      double smoothing = 0.0;
      for (int k = 0; k < dim; k++) smoothing += 2 * distance[k] * normal[k];
      if (smoothing < 0.0) {  // when nodes are ordered backwards
        normal = -normal;
        smoothing = -smoothing;
      }
      if (mesh.on_exterior_boundary(Wonton::FACE, faces[j]) and boundary_factor > 0.) {
        h[j][dim] = boundary_factor * smoothing;
      } else {
        h[j][dim] = smoothing_factor * smoothing;
      }
    }
    smoothing_lengths[i] = h;
  }

  // get extent information
  for (int i = 0; i < ncells; i++) {
    // calculate the min and max of coordinates
    std::vector<int> node_indices;
    mesh.cell_get_nodes(i, &node_indices);
    Wonton::Point<dim> cmin, cmax, first_coords;
    mesh.node_get_coordinates(node_indices[0], &first_coords);

    for (int k = 0; k < dim; k++) {
      cmin[k] = first_coords[k];
      cmax[k] = first_coords[k];
    }

    int const num_node_indices = node_indices.size();
    for (int j = 1; j < num_node_indices; j++) {
      Wonton::Point<dim> node_coords;
      mesh.node_get_coordinates(node_indices[j], &node_coords);
      for (int k = 0; k < dim; k++) {
        cmin[k] = std::min(node_coords[k], cmin[k]);
        cmax[k] = std::max(node_coords[k], cmax[k]);
      }
    }

    // subtract to get extents
    Wonton::Point<dim> dx;
    for (int k = 0; k < dim; k++) { dx[k] = cmax[k] - cmin[k]; }

    // multiply by smoothing_factor
    for (int k = 0; k < dim; k++) { dx[k] *= 2. * smoothing_factor; }

    extents[i] = dx;
  }
}


template<int dim, class Mesh>
void faceted_setup_cell(std::vector<Mesh*> meshes,
                        Portage::vector<std::vector<std::vector<double>>>& smoothing_lengths,
                        Portage::vector<Wonton::Point<dim>>& extents,
                        double smoothing_factor, double boundary_factor) {

  int total_cells = 0;
  int offset = 0;

  for (auto&& mesh : meshes)
    total_cells += mesh->num_owned_cells();

  smoothing_lengths.resize(total_cells);
  extents.resize(total_cells);

  for (auto&& mesh : meshes) {
    int ncells = mesh->num_owned_cells();
    Portage::vector<std::vector<std::vector<double>>> sl;
    Portage::vector<Wonton::Point<dim>> ex;
    faceted_setup_cell<dim, Mesh>(*mesh, sl, ex, smoothing_factor, boundary_factor);
    for (int i = 0; i < ncells; i++) {
      smoothing_lengths[offset + i] = sl[i];
      extents[offset + i] = ex[i];
    }
    offset += ncells;
  }
}


template<int dim, class Mesh, class State>
void faceted_setup_cell(const Mesh& mesh, const State& state,
                        std::string field, double tolerance,
                        Portage::vector<std::vector<std::vector<double>>>& smoothing_lengths,
                        Portage::vector<Wonton::Point<dim>>& extents,
                        double smoothing_factor, double boundary_factor) {

  faceted_setup_cell(mesh, smoothing_lengths, extents,
                     smoothing_factor, smoothing_factor);

  double *fval;
  const_cast<State&>(state).mesh_get_data(Wonton::CELL, field, &fval);

  for (int i = 0; i < mesh.num_owned_cells(); i++) {
    std::vector<int> faces, dirs, fcells;
    mesh.cell_get_faces_and_dirs(i, &faces, &dirs);
    std::vector<std::vector<double>> h = smoothing_lengths[i];
    int const num_faces = faces.size();
    for (int j = 0; j < num_faces; j++) {
      mesh.face_get_cells(faces[j], Wonton::PARALLEL_OWNED, &fcells);
      for (int k : fcells) {
        if (k == i) continue;
        if (std::fabs(fval[i] - fval[k]) > tolerance) {
          h[j][dim] *= (boundary_factor / smoothing_factor);
        }
      }
    }
    smoothing_lengths[i] = h;
  }
}

}}} // namespace Portage::Meshfree::Weight
#endif