Method initializeBoundaryMappings

src/Initializer/CellLocalMatrices.cpp:325–421

}

void seissol::initializer::initializeBoundaryMappings(const seissol::geometry::MeshReader& i_meshReader,
                                                       const EasiBoundary* easiBoundary,
                                                       LTSTree* io_ltsTree,
                                                       LTS* i_lts,
                                                       Lut* i_ltsLut) {
  std::vector<Element> const& elements = i_meshReader.getElements();
  std::vector<Vertex> const& vertices = i_meshReader.getVertices();

  const unsigned* ltsToMesh = i_ltsLut->getLtsToMeshLut(i_lts->material.mask);

  for (auto& layer : io_ltsTree->leaves(Ghost)) {
    auto* cellInformation = layer.var(i_lts->cellInformation);
    auto* boundary = layer.var(i_lts->boundaryMapping);

#ifdef _OPENMP
#pragma omp for schedule(static)
#endif
    for (unsigned cell = 0; cell < layer.getNumberOfCells(); ++cell) {
      const auto& element = elements[ltsToMesh[cell]];
      double const* coords[4];
      for (unsigned v = 0; v < 4; ++v) {
        coords[v] = vertices[ element.vertices[ v ] ].coords;
      }
      for (unsigned side = 0; side < 4; ++side) {
        if (cellInformation[cell].faceTypes[side] != FaceType::FreeSurfaceGravity
            && cellInformation[cell].faceTypes[side] != FaceType::Dirichlet
            && cellInformation[cell].faceTypes[side] != FaceType::Analytical) {
          continue;
        }
        // Compute nodal points in global coordinates for each side.
        real nodesReferenceData[nodal::tensor::nodes2D::Size];
        std::copy_n(nodal::init::nodes2D::Values,
                    nodal::tensor::nodes2D::Size,
                    nodesReferenceData);
        auto nodesReference = nodal::init::nodes2D::view::create(nodesReferenceData);
        auto nodes = boundary[cell][side].nodes;
        assert(nodes != nullptr);
        auto offset = 0;
        for (unsigned int i = 0; i < nodal::tensor::nodes2D::Shape[0]; ++i) {
          double nodeReference[2];
          nodeReference[0] = nodesReference(i,0);
          nodeReference[1] = nodesReference(i,1);
          // Compute the global coordinates for the nodal points.
          double xiEtaZeta[3], xyz[3];
          seissol::transformations::chiTau2XiEtaZeta(side,
                                                     nodeReference,
                                                     xiEtaZeta);
          seissol::transformations::tetrahedronReferenceToGlobal(coords[0],
                                                                 coords[1],
                                                                 coords[2],
                                                                 coords[3],
                                                                 xiEtaZeta,
                                                                 xyz);
          nodes[offset++] = xyz[0];
          nodes[offset++] = xyz[1];
          nodes[offset++] = xyz[2];
        }