Commit c6836617 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

Merge branch 'Add-PointInCell' into 'feature'

Add point in cell

See merge request mpp/mpp!399
parents ec2d317b 30ff1044
......@@ -71,6 +71,10 @@ public:
return (*this)->second->GlobalToLocal(z);
}
bool PointInCell(const Point &z) const {
return (*this)->second->PointInCell(z);
}
Point operator[](const Point &z) const {
return (*this)->second->LocalToGlobal(z);
}
......
......@@ -166,8 +166,13 @@ public:
virtual Point FaceLocalToGlobal(int face, const Point &) const = 0;
virtual Point LocalToGlobal(const Point &) const = 0;
virtual Point GlobalToLocal(const Point &) const = 0;
virtual bool PointInCell(const Point &) const {
Exit("Not implemented")
}
virtual Transformation GetTransformation(const Point &q = Origin) const = 0;
#ifdef BUILD_IA
......@@ -228,19 +233,10 @@ public:
return Corner(edgecorner(edge, corner));
};
/**
* Returns the index of the specified corner.
*/
virtual short edgecorner(unsigned short edge, unsigned short corner) const = 0;
/**
* Returns the amount of (d-1) dimensional faces of the cell.
*/
virtual int Faces() const = 0;
/**
* Returns the center of the (d-1) dimensional Face associated with the given integer
*/
virtual Point Face(int) const = 0;
virtual short FarFace(int) const { return -1; }
......@@ -282,18 +278,15 @@ public:
virtual double LocalFaceArea(int) const = 0;
/// Returns the area of the given face
virtual double FaceArea(int face) const = 0;
virtual const vector<Rule> &Refine(vector<Point> &) const = 0;
/// Returns the dimension of the cell
virtual int dim() const = 0;
/// Returns true, if the cell is planar
virtual bool plane() const = 0;
/// Returns the amount of children this cell generates when refined
virtual int Children() const = 0;
/// Returns the center point of the selected child
......
......@@ -139,4 +139,9 @@ Point Intval::GlobalToLocal(const Point &z) const {
return Point (x,0.0,0.0);
}
bool Intval::PointInCell(const Point &z) const {
return (abs(z[0] - Center()[0]) < 0.5 * abs(Corner(0)[0] - Corner(1)[0]) + GeometricTolerance);
}
......@@ -83,6 +83,9 @@ public:
Point Child(int i) const override;
Point GlobalToLocal(const Point &point) const override;
virtual bool PointInCell(const Point &z) const override;
};
#endif //INTVAL_H
......@@ -4,11 +4,11 @@
Quadrilateral::Quadrilateral(const vector<Point> &z,
short sd,
bool checkOrientation)
: Cell(z, sd), midpoint(0.25 * (corners[0] + corners[1] + corners[2] + corners[3])) {
: Cell(z, sd), midPoint(0.25 * (corners[0] + corners[1] + corners[2] + corners[3])) {
if (checkOrientation) this->checkOrientation(corners);
}
Point Quadrilateral::Center() const { return midpoint; }
Point Quadrilateral::Center() const { return midPoint; }
CELLTYPE Quadrilateral::Type() const { return QUADRILATERAL; }
......@@ -32,8 +32,8 @@ Point Quadrilateral::FaceLocalToGlobal(int face, const Point &local) const {
}
Point Quadrilateral::LocalToGlobal(const Point &z) const {
Point v = corners[0] + z[0] * (corners[1]-corners[0]);
Point w = corners[3] + z[0] * (corners[2]-corners[3]);
Point v = corners[0] + z[0] * (corners[1] - corners[0]);
Point w = corners[3] + z[0] * (corners[2] - corners[3]);
return v + z[1] * (w - v);
//TODO: Compare Performance
......@@ -79,7 +79,8 @@ Point Quadrilateral::Edge(int i) const {
short Quadrilateral::edgecorner(unsigned short i, unsigned short j) const {
if (i < Edges() && j < 2)
return (i + j) % 4;
else THROW("Not Implemented")
else
THROW("Not Implemented")
}
int Quadrilateral::Faces() const { return 4; }
......@@ -103,7 +104,8 @@ short Quadrilateral::FaceCorners(int i) const { return 2; }
short Quadrilateral::facecorner(unsigned short i, unsigned short j) const {
if (i < Faces() && j < FaceCorners(i))
return (i + j) % 4;
else THROW("Not Implemented")
else
THROW("Not Implemented")
}
short Quadrilateral::FaceEdges(int i) const { return 1; }
......@@ -149,7 +151,7 @@ Point Quadrilateral::Child(int i) const {
const Quadrilateral *Quadrilateral::ReferenceQuadrilateral() {
static const Quadrilateral
*refQuad = new Quadrilateral(Points(4, Quadrilateral::LocalCornersQuad()), -1);
*refQuad = new Quadrilateral(Points(4, Quadrilateral::LocalCornersQuad()), -1);
return refQuad;
}
......@@ -178,6 +180,24 @@ const Point *Quadrilateral::LocalCenterQuad() {
return localCenter;
}
bool Quadrilateral::PointInCell(const Point &z) const {
Point diff = Corner(0) - Corner(1);
double edgelength = std::max(abs(diff[0]), abs(diff[1]));
double l1distx = abs(z[0] - midPoint[0]);
double l1radius = edgelength * 0.5;
if (l1distx < l1radius + GeometricTolerance) {
double l1disty = abs(z[1] - midPoint[1]);
if (l1disty < l1radius + GeometricTolerance) {
if (abs(z[2] - midPoint[2]) < l1radius + GeometricTolerance) return true;
else return false;
} else {
return false;
}
}
return false;
}
Point Quadrilateral2::LocalToGlobal(const Point &z) const {
return ((1 - z[0]) * (1 - z[1]) * (1 - 2 * z[0] - 2 * z[1])) * corners[0]
+ (-z[0] * (1 - z[1]) * (1 - 2 * z[0] + 2 * z[1])) * corners[1]
......
......@@ -34,7 +34,7 @@ class Quadrilateral : public Cell {
}
private:
Point midpoint;
Point midPoint;
public:
Quadrilateral(const vector<Point> &z,
......@@ -75,6 +75,8 @@ public:
virtual Point LocalToGlobal(const Point &z) const override;
bool PointInCell(const Point &z) const override;
virtual Transformation GetTransformation(const Point &x) const override;
int Edges() const override;
......
#include "DGElement.hpp"
[[deprecated("Simply Use GlobalToLocal of Cell")]]
Point AcousticDGElement::GlobalToLocal(const Point &z) const {
switch(C.Type()) {
case INTERVAL: {
Point p(z);
p[0] -= C.Corner(0)[0];
double length = abs(C.Corner(0)[0] - C.Corner(1)[0]);
p /= length;
return p;
}
case TRIANGLE: {
Tensor T;
T[0][0] = C.Corner(1)[0] - C.Corner(0)[0];
T[0][1] = C.Corner(2)[0] - C.Corner(0)[0];
T[1][0] = C.Corner(1)[1] - C.Corner(0)[1];
T[1][1] = C.Corner(2)[1] - C.Corner(0)[1];
T[0][2] = 0;
T[1][2] = 0;
T[2][0] = 0;
T[2][1] = 0;
T[2][2] = 1;
Tensor IT = Invert(T);
return IT * (z - C.Corner(0));
}
case TETRAHEDRON: Exit("not implemented");
case QUADRILATERAL: {
if (inverseTensor == false) {
IT = new Tensor;
(*IT)[0][0] = C.Corner(1)[0] - C.Corner(0)[0];
(*IT)[0][1] = C.Corner(3)[0] - C.Corner(0)[0];
(*IT)[1][0] = C.Corner(1)[1] - C.Corner(0)[1];
(*IT)[1][1] = C.Corner(3)[1] - C.Corner(0)[1];
(*IT)[0][2] = 0;
(*IT)[1][2] = 0;
(*IT)[2][0] = 0;
(*IT)[2][1] = 0;
(*IT)[2][2] = 1;
IT->Invert();
inverseTensor = true;
}
return (*IT) * (z - C.Corner(0));
}
case HEXAHEDRON: Exit("not implemented");
default: Exit("not implemented");
}
Exit("Not implemented!");
return Origin;
return C.GlobalToLocal(z);
}
AcousticDGElement::AcousticDGElement(const VectorMatrixBase &base,
......
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