Loading...
Searching...
No Matches
boundarygrid.hh
Go to the documentation of this file.
//==============================================================================
//==============================================================================
#ifndef BOUNDARYGRID_HH_
#define BOUNDARYGRID_HH_
#include <opm/common/utility/platform_dependent/disable_warnings.h>
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/defaultgridview.hh>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <vector>
namespace Opm {
namespace Elasticity {
class BoundaryGrid {
public:
typedef Dune::FieldVector<double,3> GlobalCoordinate;
typedef Dune::FieldVector<double,2> FaceCoord;
// in natural order along the first direction
static BoundaryGrid uniform(const FaceCoord& min, const FaceCoord& max,
int k1, int k2, bool dc=false);
BoundaryGrid() {}
virtual ~BoundaryGrid() {}
// on a boundary
class Quad;
class Vertex {
public:
Vertex() : i(-1), c(0), q(0), fixed(false) {}
int i;
FaceCoord c;
Quad* q;
bool fixed;
bool operator==(const Vertex& v2)
{
return hypot(v2.c[0]-c[0],v2.c[1]-c[1]) < 1.e-8;
}
};
class Quad {
public:
Quad()
{
v[0].i = v[1].i = v[2].i = v[3].i = 0;
v[0].c = v[1].c = v[2].c = v[3].c = 0.f;
}
FaceCoord pos(double xi, double eta) const;
std::vector<double> evalBasis(double xi, double eta) const;
Vertex v[4];
FaceCoord bb[2];
protected:
friend std::ostream& operator <<(std::ostream& os, const Quad& q)
{
os << "Nodes: " << q.v[0].i << "," << q.v[1].i << "," << q.v[2].i << "," << q.v[3].i << std::endl;
os << "(" << q.v[0].c << ")(" << q.v[1].c << ")(" << q.v[2].c << ")(" << q.v[3].c << ")";
return os;
}
};
void add(const Quad& quad);
void addToColumn(size_t col, const Quad& quad)
{
if (col >= colGrids.size())
colGrids.resize(col+1);
colGrids[col].push_back(quad);
}
Quad& operator[](int index)
{
return grid[index];
}
const Quad& operator[](int index) const
{
return grid[index];
}
const Quad& getQuad(int col, int index) const
{
return colGrids[col][index];
}
size_t size() const
{
return grid.size();
}
size_t colSize(int i) const
{
return colGrids[i].size();
}
size_t totalNodes() const
{
return nodes;
}
bool isFixed(int node) const
{
return fixNodes[node];
}
bool find(Vertex& res, const Vertex& coord) const;
static void extract(FaceCoord& res,
const GlobalCoordinate& coord, int dir);
static void extract(Vertex& res,
const Dune::FieldVector<double,3>& coord, int dir);
struct VertexLess {
explicit VertexLess(int comp) : dir(comp) {}
bool operator()(const Vertex& q1, const Vertex& q2)
{
if (dir >= 0)
return q1.c[dir] < q2.c[dir];
return q1.i < q2.i;
}
int dir;
};
struct BoundedPredicate {
explicit BoundedPredicate(const FaceCoord& coord_) : coord(coord_) {}
bool operator()(const Quad& q)
{
double eps = 1.e-8;
return (coord[0] >= q.bb[0][0]-eps &&
coord[0] <= q.bb[1][0]+eps &&
coord[1] >= q.bb[0][1]-eps &&
coord[1] <= q.bb[1][1]+eps);
}
FaceCoord coord;
};
protected:
std::vector<Quad> grid;
std::vector<std::vector<Quad> > colGrids;
std::vector<bool> fixNodes;
size_t nodes;
friend std::ostream& operator <<(std::ostream& os, const BoundaryGrid& g)
{
for (size_t i=0;i<g.size();++i)
os << g[i] << std::endl;
return os;
}
bool bilinearSolve(double epsilon, double order,
const double* A, const double* B,
std::vector<double>& X,
std::vector<double>& Y) const;
bool cubicSolve(double eps, double A, double B, double C,
double D, std::vector<double>& X) const;
int Q4inv(FaceCoord& res, const Quad& q, const FaceCoord& coord,
double epsZero, double epsOut) const;
};
template<int mydim, int cdim, class GridImp>
class HexGeometry
{
};
template<int cdim, class GridImp>
class HexGeometry<2, cdim, GridImp>
{
public:
enum { dimension = 2};
enum { mydimension = 2};
enum { coorddimension = cdim };
enum {dimensionworld = 2 };
typedef double ctype;
typedef Dune::FieldVector<ctype,mydimension> LocalCoordinate;
typedef Dune::FieldVector<ctype,coorddimension> GlobalCoordinate;
typedef Dune::FieldMatrix<ctype,coorddimension,mydimension> Jacobian;
typedef Dune::FieldMatrix<ctype,mydimension,coorddimension> JacobianTransposed;
HexGeometry(const BoundaryGrid::Quad& q, const GridImp& gv, int dir)
{
using LeafGridView = Dune::GridView<Dune::DefaultLeafGridViewTraits<GridImp>>;
Dune::MultipleCodimMultipleGeomTypeMapper<LeafGridView> mapper(gv.leafGridView(), Dune::mcmgVertexLayout());
typename LeafGridView::template Codim<3>::Iterator start=gv.leafGridView().template begin<3>();
const typename LeafGridView::template Codim<3>::Iterator itend = gv.leafGridView().template end<3>();
for (int i=0;i<4;++i) {
typename LeafGridView::template Codim<3>::Iterator it=start;
for (; it != itend; ++it) {
if (mapper.index(*it) == q.v[i].i)
break;
}
BoundaryGrid::extract(c[i],it->geometry().corner(0),dir);
}
m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
}
explicit HexGeometry(const BoundaryGrid::Quad& q)
{
for (int i=0;i<4;++i)
c[i] = q.v[i].c;
m_volume = (c[1][0]-c[0][0])*(c[2][1]-c[0][1]);
}
Dune::GeometryType type() const
{
Dune::GeometryType t;
t = Dune::GeometryTypes::cube(mydimension);
return t;
}
int corners() const
{
return 4;
}
ctype volume() const
{
return m_volume;
}
GlobalCoordinate center() const
{
LocalCoordinate local;
local = .5f;
return Global(local);
}
GlobalCoordinate corner(int cor) const
{
return c[cor];
}
GlobalCoordinate global(const LocalCoordinate& local) const
{
// uvw = { (1-u, 1-v, 1-w), (u, v, w) }
LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
uvw[0] -= local;
// Access pattern for uvw matching ordering of corners.
const int pat[4][2] = {{ 0, 0 },
{ 1, 0 },
{ 0, 1 },
{ 1, 1 }};
GlobalCoordinate xyz(0.0);
for (int i = 0; i < 4; ++i) {
GlobalCoordinate corner_contrib = corner(i);
double factor = 1.0;
for (int j = 0; j < 2; ++j) {
factor *= uvw[pat[i][j]][j];
}
corner_contrib *= factor;
xyz += corner_contrib;
}
return xyz;
}
LocalCoordinate local(const GlobalCoordinate& y) const
{
const ctype epsilon = 1e-10;
auto refElement = Dune::ReferenceElements<ctype,2>::cube();
LocalCoordinate x = refElement.position(0,0);
LocalCoordinate dx;
do {
// DF^n dx^n = F^n, x^{n+1} -= dx^n
JacobianTransposed JT = jacobianTransposed(x);
GlobalCoordinate z = global(x);
z -= y;
Dune::Impl::FieldMatrixHelper<double>::template xTRightInvA<2, 2>(JT, z, dx );
x -= dx;
} while (dx.two_norm2() > epsilon*epsilon);
return x;
}
const Dune::FieldMatrix<ctype, mydimension, coorddimension>
jacobianTransposed(const LocalCoordinate& local) const
{
// uvw = { (1-u, 1-v, (u, v) }
LocalCoordinate uvw[2] = { LocalCoordinate(1.0), local };
uvw[0] -= local;
// Access pattern for uvw matching ordering of corners.
const int pat[4][3] = {{ 0, 0 },
{ 1, 0 },
{ 0, 1 },
{ 1, 1 }};
Dune::FieldMatrix<ctype, mydimension, coorddimension> Jt(0.0);
for (int i = 0; i < 4; ++i) {
for (int deriv = 0; deriv < 2; ++deriv) {
// This part contributing to dg/du_{deriv}
double factor = 1.0;
for (int j = 0; j < 2; ++j) {
factor *= (j != deriv) ? uvw[pat[i][j]][j]
: (pat[i][j] == 0 ? -1.0 : 1.0);
}
GlobalCoordinate corner_contrib = corner(i);
corner_contrib *= factor;
Jt[deriv] += corner_contrib; // using FieldMatrix row access.
}
}
return Jt;
}
const Dune::FieldMatrix<ctype, coorddimension, mydimension>
jacobianInverseTransposed(const LocalCoordinate& local) const
{
Dune::FieldMatrix<ctype, coorddimension, mydimension> Jti = jacobianTransposed(local);
Jti.invert();
return Jti;
}
ctype integrationElement(const LocalCoordinate& local) const
{
Dune::FieldMatrix<ctype, coorddimension, mydimension> Jt = jacobianTransposed(local);
return Dune::Impl::FieldMatrixHelper<double>::template sqrtDetAAT<2, 2>(Jt);
}
private:
GlobalCoordinate c[4];
ctype m_volume;
};
BoundaryGrid::Vertex minXminY(const std::vector<BoundaryGrid::Vertex>& in);
BoundaryGrid::Vertex maxXminY(const std::vector<BoundaryGrid::Vertex>& in);
BoundaryGrid::Vertex maxXmaxY(const std::vector<BoundaryGrid::Vertex>& in);
BoundaryGrid::Vertex minXmaxY(const std::vector<BoundaryGrid::Vertex>& in);
}
}
#endif
Opm::Elasticity::BoundaryGrid::Quad
A class describing a linear, quadrilateral element.
Definition boundarygrid.hh:86
Opm::Elasticity::BoundaryGrid::Quad::Quad
Quad()
Default constructor.
Definition boundarygrid.hh:89
Opm::Elasticity::BoundaryGrid::Quad::v
Vertex v[4]
Vertices.
Definition boundarygrid.hh:106
Opm::Elasticity::BoundaryGrid::Quad::bb
FaceCoord bb[2]
Bounding box.
Definition boundarygrid.hh:108
Opm::Elasticity::BoundaryGrid::Quad::operator<<
friend std::ostream & operator<<(std::ostream &os, const Quad &q)
Print to a stream.
Definition boundarygrid.hh:111
Opm::Elasticity::BoundaryGrid::Quad::evalBasis
std::vector< double > evalBasis(double xi, double eta) const
Evaluate the basis functions.
Definition boundarygrid.cpp:340
Opm::Elasticity::BoundaryGrid::Quad::pos
FaceCoord pos(double xi, double eta) const
Return the physical coordinates corresponding to the given local coordinates.
Definition boundarygrid.cpp:329
Opm::Elasticity::BoundaryGrid::Vertex
A class describing a 2D vertex.
Definition boundarygrid.hh:62
Opm::Elasticity::BoundaryGrid::Vertex::operator==
bool operator==(const Vertex &v2)
Comparison operator.
Definition boundarygrid.hh:80
Opm::Elasticity::BoundaryGrid::Vertex::q
Quad * q
The quad containing the vertex (if search has been done).
Definition boundarygrid.hh:74
Opm::Elasticity::BoundaryGrid::Vertex::fixed
bool fixed
Whether or not this node is fixed.
Definition boundarygrid.hh:77
Opm::Elasticity::BoundaryGrid::Vertex::c
FaceCoord c
Coordinates of the vertex (2D).
Definition boundarygrid.hh:71
Opm::Elasticity::BoundaryGrid::Vertex::i
int i
Index of the vertex.
Definition boundarygrid.hh:68
Opm::Elasticity::BoundaryGrid::Vertex::Vertex
Vertex()
Default constructor.
Definition boundarygrid.hh:65
Opm::Elasticity::BoundaryGrid::grid
std::vector< Quad > grid
Our quadrilateral elements.
Definition boundarygrid.hh:235
Opm::Elasticity::BoundaryGrid::fixNodes
std::vector< bool > fixNodes
Whether or not a given node is marked as fixed.
Definition boundarygrid.hh:239
Opm::Elasticity::BoundaryGrid::cubicSolve
bool cubicSolve(double eps, double A, double B, double C, double D, std::vector< double > &X) const
Solves the cubic equation A*x^3+B*x^2+C*x+D=0.
Definition boundarygrid.cpp:270
Opm::Elasticity::BoundaryGrid::operator[]
Quad & operator[](int index)
Obtain a reference to a quad.
Definition boundarygrid.hh:133
Opm::Elasticity::BoundaryGrid::~BoundaryGrid
virtual ~BoundaryGrid()
Default (empty) destructor.
Definition boundarygrid.hh:55
Opm::Elasticity::BoundaryGrid::GlobalCoordinate
Dune::FieldVector< double, 3 > GlobalCoordinate
A coordinate on the underlying 3D grid.
Definition boundarygrid.hh:36
Opm::Elasticity::BoundaryGrid::totalNodes
size_t totalNodes() const
Return the total number of nodes on the grid when known.
Definition boundarygrid.hh:164
Opm::Elasticity::BoundaryGrid::nodes
size_t nodes
Total number of nodes on grid.
Definition boundarygrid.hh:242
Opm::Elasticity::BoundaryGrid::BoundaryGrid
BoundaryGrid()
Default (empty) constructor.
Definition boundarygrid.hh:52
Opm::Elasticity::BoundaryGrid::operator[]
const Quad & operator[](int index) const
Obtain a const reference to a quad.
Definition boundarygrid.hh:141
Opm::Elasticity::BoundaryGrid::operator<<
friend std::ostream & operator<<(std::ostream &os, const BoundaryGrid &g)
Print to a stream.
Definition boundarygrid.hh:245
Opm::Elasticity::BoundaryGrid::add
void add(const Quad &quad)
Add a quad to the grid.
Definition boundarygrid.cpp:85
Opm::Elasticity::BoundaryGrid::find
bool find(Vertex &res, const Vertex &coord) const
Locate the position of a vertex on the grid.
Definition boundarygrid.cpp:117
Opm::Elasticity::BoundaryGrid::extract
static void extract(FaceCoord &res, const GlobalCoordinate &coord, int dir)
Helper function for extracting given 2D coordinates from a 3D coordinate.
Definition boundarygrid.cpp:70
Opm::Elasticity::BoundaryGrid::bilinearSolve
bool bilinearSolve(double epsilon, double order, const double *A, const double *B, std::vector< double > &X, std::vector< double > &Y) const
Solves a bi-linear set of equations in x and y.
Definition boundarygrid.cpp:213
Opm::Elasticity::BoundaryGrid::FaceCoord
Dune::FieldVector< double, 2 > FaceCoord
A coordinate on the quad grid.
Definition boundarygrid.hh:38
Opm::Elasticity::BoundaryGrid::Q4inv
int Q4inv(FaceCoord &res, const Quad &q, const FaceCoord &coord, double epsZero, double epsOut) const
Find the local coordinates of a given point within a given quad.
Definition boundarygrid.cpp:150
Opm::Elasticity::BoundaryGrid::isFixed
bool isFixed(int node) const
Check if a given node is marked as fixed.
Definition boundarygrid.hh:172
Opm::Elasticity::BoundaryGrid::uniform
static BoundaryGrid uniform(const FaceCoord &min, const FaceCoord &max, int k1, int k2, bool dc=false)
Establish an uniform quad grid.
Definition boundarygrid.cpp:24
Opm::Elasticity::BoundaryGrid::size
size_t size() const
Obtain the number of quads in the grid.
Definition boundarygrid.hh:152
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::Jacobian
Dune::FieldMatrix< ctype, coorddimension, mydimension > Jacobian
Type of Jacobian matrix.
Definition boundarygrid.hh:319
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::HexGeometry
HexGeometry(const BoundaryGrid::Quad &q, const GridImp &gv, int dir)
Construct integration element extracted from a 3D grid.
Definition boundarygrid.hh:328
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::corners
int corners() const
Returns number of corners.
Definition boundarygrid.hh:364
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::center
GlobalCoordinate center() const
Returns center of quadrilateral.
Definition boundarygrid.hh:376
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::local
LocalCoordinate local(const GlobalCoordinate &y) const
Map from global coordinates to local coordinates.
Definition boundarygrid.hh:417
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::jacobianInverseTransposed
const Dune::FieldMatrix< ctype, coorddimension, mydimension > jacobianInverseTransposed(const LocalCoordinate &local) const
Returns the inverse, transposed Jacobian.
Definition boundarygrid.hh:467
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::LocalCoordinate
Dune::FieldVector< ctype, mydimension > LocalCoordinate
Domain type.
Definition boundarygrid.hh:313
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::type
Dune::GeometryType type() const
Returns entity type (a 2D cube).
Definition boundarygrid.hh:356
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::volume
ctype volume() const
Returns volume (area) of quadrilateral.
Definition boundarygrid.hh:370
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::ctype
double ctype
Coordinate element type.
Definition boundarygrid.hh:310
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::GlobalCoordinate
Dune::FieldVector< ctype, coorddimension > GlobalCoordinate
Range type.
Definition boundarygrid.hh:316
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::JacobianTransposed
Dune::FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Type of transposed Jacobian matrix.
Definition boundarygrid.hh:322
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::corner
GlobalCoordinate corner(int cor) const
Returns coordinates to requested corner.
Definition boundarygrid.hh:385
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::integrationElement
ctype integrationElement(const LocalCoordinate &local) const
Returns the integration element (|J'*J|)^(1/2).
Definition boundarygrid.hh:476
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::jacobianTransposed
const Dune::FieldMatrix< ctype, mydimension, coorddimension > jacobianTransposed(const LocalCoordinate &local) const
Return the transposed jacobian.
Definition boundarygrid.hh:437
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::HexGeometry
HexGeometry(const BoundaryGrid::Quad &q)
Construct integration element.
Definition boundarygrid.hh:348
Opm::Elasticity::HexGeometry< 2, cdim, GridImp >::global
GlobalCoordinate global(const LocalCoordinate &local) const
Map from local coordinates to global coordinates.
Definition boundarygrid.hh:392
Opm::Elasticity::HexGeometry
Geometry class for general hexagons.
Definition boundarygrid.hh:289
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43
Opm::Elasticity::BoundaryGrid::BoundedPredicate::operator()
bool operator()(const Quad &q)
The comparison operator.
Definition boundarygrid.hh:221
Opm::Elasticity::BoundaryGrid::BoundedPredicate::coord
FaceCoord coord
The coordinates to check.
Definition boundarygrid.hh:231
Opm::Elasticity::BoundaryGrid::BoundedPredicate::BoundedPredicate
BoundedPredicate(const FaceCoord &coord_)
Default constructor.
Definition boundarygrid.hh:218
Opm::Elasticity::BoundaryGrid::VertexLess::VertexLess
VertexLess(int comp)
Default constructor.
Definition boundarygrid.hh:200
Opm::Elasticity::BoundaryGrid::VertexLess::dir
int dir
Compare using this direction, if -1 compare using index values.
Definition boundarygrid.hh:211
Opm::Elasticity::BoundaryGrid::VertexLess::operator()
bool operator()(const Vertex &q1, const Vertex &q2)
The comparison operator.
Definition boundarygrid.hh:203