Commit 9d9e7984 authored by Timo Koch's avatar Timo Koch

[!57] Fix/geometryinoutside 1d grids

Merge branch 'fix/geometryinoutside-1d-grids' into 'master'

See merge request [!57]

  [!57]: Noneextensions/dune-foamgrid/merge_requests/57
parents 26dd3047 f1805568
Pipeline #18174 passed with stage
in 2 minutes and 30 seconds
......@@ -117,6 +117,23 @@ namespace Dune {
return InteriorEntity;
}
/** \brief Compute local cordinates from global ones.
* \param coord The global coordinates.
* \return The corresponding local coordinates within the element.
*/
FieldVector<double, 1> globalToLocal(const FieldVector<double, dimworld>& coord) const
{
const auto diff = vertex_[1]->pos_ - vertex_[0]->pos_;
const double eps = diff.two_norm()*std::numeric_limits<double>::epsilon();
for (std::size_t dimIdx = 0; dimIdx < dimworld; ++dimIdx)
if (std::abs(diff[dimIdx]) > eps)
return (coord[dimIdx] - vertex_[0]->pos_[dimIdx]) / diff[dimIdx];
DUNE_THROW(Dune::GridError, "Global to local mapping failed because element is degenerated.");
}
/** \brief Return level index of sub entity with codim = cc and local number i
*/
int subLevelIndex (int i, unsigned int codim) const {
......
......@@ -477,27 +477,20 @@ public:
//! intersection of codimension 1 of this neighbor with element where iteration started.
//! Here returned element is in LOCAL coordinates of neighbor
//! In the LevelIntersection we know that the intersection is conforming
LocalGeometry geometryInOutside (std::size_t neighborIndex = 0) const {
// Get two vertices of the intersection
const auto refElement = ReferenceElements<double, dimgrid>::general(this->center_->type());
std::array<FoamGridEntityImp<0, dimgrid, dimworld>*, dimgrid> vtx;
for (std::size_t idx = 0; idx < dimgrid; ++idx)
vtx[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)];
//! In the LeafIntersection the intersection might be non-conforming
//! For surface grids with local non-conforming adaptivity the geometryInOutside is undefined
LocalGeometry geometryInOutside (std::size_t neighborIndex = 0) const
{
// Get reference element
const auto refElement = Dune::ReferenceElements<double, dimgrid>::general(this->center_->type());
// Map the global position of the intersection vertices to the local position in the neighbor
std::vector<FieldVector<double, dimgrid> > coordinates(dimgrid);
// Find the intersection vertices in local numbering of the outside element
// That way we get the local orientation correctly.
const auto refElementOther = ReferenceElements<double, dimgrid>::general((*this->neighbor_)->type());
for (std::size_t j=0; j< dimgrid; j++)
for (int i=0; i<refElementOther.size(dimgrid); i++)
if (vtx[j] == (*this->neighbor_)->vertex_[refElementOther.subEntity(0, 0, i, dimgrid)])
coordinates[j] = refElement.position(refElement.subEntity(0, 0, i, dimgrid), dimgrid);
for (std::size_t vIdx = 0; vIdx < dimgrid; ++vIdx)
{
const auto vIdxInElement = refElement.subEntity(this->facetIndex_, 1, vIdx, dimgrid);
coordinates[vIdx] = (*this->neighbor_)->globalToLocal(this->center_->vertex_[vIdxInElement]->pos_);
}
geometryInOutside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment