Loading...
Searching...
No Matches
codim1extractor.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-GPL-2.0-only-with-dune-grid-glue-exception
/*
* Filename: codim1extractor.hh
* Version: 1.0
* Created on: Jun 23, 2009
* Author: Oliver Sander, Christian Engwer
* ---------------------------------
* Project: dune-grid-glue
* Description: class for grid extractors extracting surface grids
*
*/
#ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
#define DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
#include "extractor.hh"
#include <array>
#include <deque>
#include <functional>
#include <dune/common/deprecated.hh>
#include <dune/common/version.hh>
#include <dune/grid-glue/common/crossproduct.hh>
namespace Dune {
namespace GridGlue {
template<typename GV>
class Codim1Extractor : public Extractor<GV,1>
{
public:
/* E X P O R T E D T Y P E S A N D C O N S T A N T S */
using Extractor<GV,1>::dimworld;
using Extractor<GV,1>::dim;
using Extractor<GV,1>::codim;
using Extractor<GV,1>::cube_corners;
typedef typename Extractor<GV,1>::IndexType IndexType;
static constexpr int simplex_corners = dim;
typedef GV GridView;
typedef typename GV::Grid::ctype ctype;
typedef Dune::FieldVector<ctype, dimworld> Coords;
typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
typedef typename GV::Traits::template Codim<0>::Entity Element;
typedef std::function<bool(const Element&, unsigned int subentity)> Predicate;
// import typedefs from base class
typedef typename Extractor<GV,1>::SubEntityInfo SubEntityInfo;
typedef typename Extractor<GV,1>::ElementInfo ElementInfo;
typedef typename Extractor<GV,1>::VertexInfo VertexInfo;
typedef typename Extractor<GV,1>::CoordinateInfo CoordinateInfo;
typedef typename Extractor<GV,1>::VertexInfoMap VertexInfoMap;
public:
/* C O N S T R U C T O R S A N D D E S T R U C T O R S */
Codim1Extractor(const GV& gv, const Predicate& predicate)
: Extractor<GV,1>(gv)
{
std::cout << "This is Codim1Extractor on a <" << dim
<< "," << dimworld << "> grid!"
<< std::endl;
update(predicate);
}
private:
void update(const Predicate& predicate);
};
template<typename GV>
void Codim1Extractor<GV>::update(const Predicate& predicate)
{
// free everything there is in this object
this->clear();
// In this first pass iterate over all entities of codim 0.
// For each codim 1 intersection check if it is part of the boundary and if so,
// get its corner vertices, find resp. store them together with their associated index,
// and remember the indices of the boundary faces' corners.
{
// several counter for consecutive indexing are needed
int simplex_index = 0;
int vertex_index = 0;
IndexType eindex = 0; // suppress warning
// needed later for insertion into a std::set which only
// works with const references
// a temporary container where newly acquired face
// information can be stored at first
std::deque<SubEntityInfo> temp_faces;
// iterate over interior codim 0 elements on the grid
for (const auto& elmt : elements(this->gv_, Partitions::interior))
{
Dune::GeometryType gt = elmt.type();
// if some face is part of the surface add it!
if (elmt.hasBoundaryIntersections())
{
// add an entry to the element info map, the index will be set properly later,
// whereas the number of faces is already known
eindex = this->cellMapper_.index(elmt);
this->elmtInfo_.emplace(eindex, ElementInfo(simplex_index, elmt, 0));
// now add the faces in ascending order of their indices
// (we are only talking about 1-4 faces here, so O(n^2) is ok!)
for (const auto& in : intersections(this->gv_, elmt))
{
// Stop only at selected boundary faces
if (!in.boundary() or !predicate(elmt, in.indexInInside()))
continue;
const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
// get the corner count of this face
const int face_corners = refElement.size(in.indexInInside(), 1, dim);
// now we only have to care about the 3D case, i.e. a triangle face can be
// inserted directly whereas a quadrilateral face has to be divided into two triangles
switch (face_corners)
{
case 2 :
case 3:
{
// we have a simplex here
// register the additional face(s)
this->elmtInfo_.at(eindex).faces++;
// add a new face to the temporary collection
temp_faces.emplace_back(eindex, in.indexInInside(),
Dune::GeometryTypes::simplex(dim-codim));
std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
// try for each of the faces vertices whether it is already inserted or not
for (int i = 0; i < face_corners; ++i)
{
// get the number of the vertex in the parent element
int vertex_number = refElement.subEntity(in.indexInInside(), 1, i, dim);
// get the vertex pointer and the index from the index set
const Vertex vertex = elmt.template subEntity<dim>(vertex_number);
cornerCoords[i] = vertex.geometry().corner(0);
IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
// remember the vertex' number in parent element's vertices
temp_faces.back().corners[i].num = vertex_number;
// if the vertex is not yet inserted in the vertex info map
// it is a new one -> it will be inserted now!
typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
if (vimit == this->vtxInfo_.end())
{
// insert into the map
this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
// remember the vertex as a corner of the current face in temp_faces
temp_faces.back().corners[i].idx = vertex_index;
// increase the current index
vertex_index++;
}
else
{
// only insert the index into the simplices array
temp_faces.back().corners[i].idx = vimit->second.idx;
}
}
// Now we have the correct vertices in the last entries of temp_faces, but they may
// have the wrong orientation. We want them to be oriented such that all boundary edges
// point in the counterclockwise direction. Therefore, we check the orientation of the
// new face and possibly switch the two vertices.
FieldVector<ctype,dimworld> realNormal = in.centerUnitOuterNormal();
// Compute segment normal
FieldVector<ctype,dimworld> reconstructedNormal;
if (dim==2) // boundary face is a line segment
{
reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
} else { // boundary face is a triangle
FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
reconstructedNormal = crossProduct(segment1, segment2);
}
reconstructedNormal /= reconstructedNormal.two_norm();
if (realNormal * reconstructedNormal < 0.0)
std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
// now increase the current face index
simplex_index++;
break;
}
case 4 :
{
assert(dim == 3 && cube_corners == 4);
// we have a quadrilateral here
std::array<unsigned int, 4> vertex_indices;
std::array<unsigned int, 4> vertex_numbers;
// register the additional face(s) (2 simplices)
this->elmtInfo_.at(eindex).faces += 2;
std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
// get the vertex pointers for the quadrilateral's corner vertices
// and try for each of them whether it is already inserted or not
for (int i = 0; i < cube_corners; ++i)
{
// get the number of the vertex in the parent element
vertex_numbers[i] = refElement.subEntity(in.indexInInside(), 1, i, dim);
// get the vertex pointer and the index from the index set
const Vertex vertex = elmt.template subEntity<dim>(vertex_numbers[i]);
cornerCoords[i] = vertex.geometry().corner(0);
IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
// if the vertex is not yet inserted in the vertex info map
// it is a new one -> it will be inserted now!
typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
if (vimit == this->vtxInfo_.end())
{
// insert into the map
this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
// remember this vertex' index
vertex_indices[i] = vertex_index;
// increase the current index
vertex_index++;
}
else
{
// only remember the vertex' index
vertex_indices[i] = vimit->second.idx;
}
}
// now introduce the two triangles subdividing the quadrilateral
// ATTENTION: the order of vertices given by "orientedSubface" corresponds to the order
// of a Dune quadrilateral, i.e. the triangles are given by 0 1 2 and 3 2 1
// add a new face to the temporary collection for the first tri
temp_faces.emplace_back(eindex, in.indexInInside(),
Dune::GeometryTypes::simplex(dim-codim));
temp_faces.back().corners[0].idx = vertex_indices[0];
temp_faces.back().corners[1].idx = vertex_indices[1];
temp_faces.back().corners[2].idx = vertex_indices[2];
// remember the vertices' numbers in parent element's vertices
temp_faces.back().corners[0].num = vertex_numbers[0];
temp_faces.back().corners[1].num = vertex_numbers[1];
temp_faces.back().corners[2].num = vertex_numbers[2];
// Now we have the correct vertices in the last entries of temp_faces, but they may
// have the wrong orientation. We want the triangle vertices on counterclockwise order,
// when viewed from the outside of the grid. Therefore, we check the orientation of the
// new face and possibly switch two vertices.
FieldVector<ctype,dimworld> realNormal = in.centerUnitOuterNormal();
// Compute segment normal
FieldVector<ctype,dimworld> reconstructedNormal = crossProduct(cornerCoords[1] - cornerCoords[0],
cornerCoords[2] - cornerCoords[0]);
reconstructedNormal /= reconstructedNormal.two_norm();
if (realNormal * reconstructedNormal < 0.0)
std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
// add a new face to the temporary collection for the second tri
temp_faces.emplace_back(eindex, in.indexInInside(),
Dune::GeometryTypes::simplex(dim-codim));
temp_faces.back().corners[0].idx = vertex_indices[3];
temp_faces.back().corners[1].idx = vertex_indices[2];
temp_faces.back().corners[2].idx = vertex_indices[1];
// remember the vertices' numbers in parent element's vertices
temp_faces.back().corners[0].num = vertex_numbers[3];
temp_faces.back().corners[1].num = vertex_numbers[2];
temp_faces.back().corners[2].num = vertex_numbers[1];
// Now we have the correct vertices in the last entries of temp_faces, but they may
// have the wrong orientation. We want the triangle vertices on counterclockwise order,
// when viewed from the outside of the grid. Therefore, we check the orientation of the
// new face and possibly switch two vertices.
// Compute segment normal
reconstructedNormal = crossProduct(cornerCoords[2] - cornerCoords[3],
cornerCoords[1] - cornerCoords[3]);
reconstructedNormal /= reconstructedNormal.two_norm();
if (realNormal * reconstructedNormal < 0.0)
std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
simplex_index+=2;
break;
}
default :
DUNE_THROW(Dune::NotImplemented, "the extractor does only work for triangle and quadrilateral faces (" << face_corners << " corners)");
break;
}
} // end loop over found surface parts
}
} // end loop over elements
std::cout << "added " << simplex_index << " subfaces\n";
// allocate the array for the face specific information...
this->subEntities_.resize(simplex_index);
// ...and fill in the data from the temporary containers
copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
}
// now first write the array with the coordinates...
this->coords_.resize(this->vtxInfo_.size());
for (const auto& vinfo : this->vtxInfo_)
{
// get a pointer to the associated info object
CoordinateInfo* current = &this->coords_[vinfo.second.idx];
// store this coordinates index // NEEDED?
current->index = vinfo.second.idx;
// store the vertex' index for the index2vertex mapping
current->vtxindex = vinfo.first;
// store the vertex' coordinates under the associated index
// in the coordinates array
const auto vtx = this->grid().entity(vinfo.second.p);
current->coord = vtx.geometry().corner(0);
}
}
} // namespace GridGlue
} // namespace Dune
#endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
crossproduct.hh
extractor.hh
extractor base class
Dune
Definition gridglue.hh:37
Dune::GridGlue
Definition gridglue.hh:38
Dune::GridGlue::intersections
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Dune::GridGlue::crossProduct
static Dune::FieldVector< T, dim > crossProduct(const Dune::FieldVector< T, dim > &a, const Dune::FieldVector< T, dim > &b)
compute cross product
Definition crossproduct.hh:15
Dune::GridGlue::Codim1Extractor::IndexType
Extractor< GV, 1 >::IndexType IndexType
Definition codim1extractor.hh:51
Dune::GridGlue::Codim1Extractor::GridView
GV GridView
Definition codim1extractor.hh:56
Dune::GridGlue::Codim1Extractor::Element
GV::Traits::template Codim< 0 >::Entity Element
Definition codim1extractor.hh:62
Dune::GridGlue::Codim1Extractor::Coords
Dune::FieldVector< ctype, dimworld > Coords
Definition codim1extractor.hh:59
Dune::GridGlue::Codim1Extractor::ctype
GV::Grid::ctype ctype
Definition codim1extractor.hh:58
Dune::GridGlue::Codim1Extractor::Codim1Extractor
Codim1Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition codim1extractor.hh:81
Dune::GridGlue::Codim1Extractor::simplex_corners
static constexpr int simplex_corners
compile time number of corners of surface simplices
Definition codim1extractor.hh:54
Dune::GridGlue::Codim1Extractor::VertexInfo
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition codim1extractor.hh:68
Dune::GridGlue::Codim1Extractor::CoordinateInfo
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition codim1extractor.hh:69
Dune::GridGlue::Codim1Extractor::ElementInfo
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition codim1extractor.hh:67
Dune::GridGlue::Codim1Extractor::SubEntityInfo
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition codim1extractor.hh:66
Dune::GridGlue::Codim1Extractor::VertexInfoMap
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition codim1extractor.hh:70
Dune::GridGlue::Codim1Extractor::Predicate
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition codim1extractor.hh:63
Dune::GridGlue::Codim1Extractor::Vertex
GV::Traits::template Codim< dim >::Entity Vertex
Definition codim1extractor.hh:61
Dune::GridGlue::IndexType
Dune::GridGlue::Extractor< GV, 1 >::dimworld
static constexpr auto dimworld
Definition extractor.hh:50
Dune::GridGlue::Extractor::IndexType
int IndexType
Definition extractor.hh:73
Dune::GridGlue::Extractor< GV, 1 >::cube_corners
static constexpr int cube_corners
Definition extractor.hh:54
Dune::GridGlue::Extractor< GV, 1 >::codim
static constexpr auto codim
Definition extractor.hh:52
Dune::GridGlue::Extractor::VertexInfoMap
std::map< IndexType, VertexInfo > VertexInfoMap
Definition extractor.hh:186
Dune::GridGlue::Extractor< GV, 1 >::Extractor
Extractor(const GV &gv)
Definition extractor.hh:225
Dune::GridGlue::Extractor< GV, 1 >::dim
static constexpr auto dim
Definition extractor.hh:51
Dune::GridGlue::Extractor::CoordinateInfo
Definition extractor.hh:94
Dune::GridGlue::Extractor::VertexInfo
simple struct holding a vertex pointer and an index
Definition extractor.hh:116
Dune::GridGlue::Extractor::ElementInfo
simple struct holding an element seed and an index
Definition extractor.hh:128
Dune::GridGlue::Extractor::SubEntityInfo
Holds some information about an element's subentity involved in a coupling.
Definition extractor.hh:147