Commit 661796b4 authored by Nicolas CHEVAUGEON's avatar Nicolas CHEVAUGEON

[xMapping] adding some virtual member to xMapping.

xMapping now include xTensor/src/xBoundingBox.

The following virtual member function of xMapping have been added.
They all have default implementation in terms of the one existing before.
They mostly simplify code when using them, and do not need to be implemented i
in derived class (but could for optimization purpose)

 virtual xtensor::xPoint eval(const xtensor::xPoint &uvw) const;
 virtual tensor deval(const xtensor::xPoint &uvw) const;
 virtual bool invert(const xtensor::xPoint &xyz, xtensor::xPoint &uvw) const;
 virtual double jacInverse(const xtensor::xPoint &uvw, tensor &invjac) const;
 virtual double detJac(const xtensor::xPoint &uvw) const;
 virtual bool inReferenceElement(const xtensor::xPoint &uvw) const;
 virtual xtensor::xPoint COG() const;
 virtual double pushBack(const xtensor::xPoint &uvw, size_t vsize, vector3d *vec) const;
 virtual double pushBack(const xtensor::xPoint &uvw, std::vector<vector3d> &vec) const;
 virtual double pushBack(const xtensor::xPoint &uvw, size_t vsize, tensor *tens) const;
 virtual xtensor::xBoundingBox boundingBox() const;

interface change in xMapping :
    virtual bool interiorCheck(const xtensor::xPoint &p, double u, double v, double w) const; changed to virtual bool interiorCheck(const xtensor::xPoint &p, xtensor::xPoint &uvw) const;

 Modified files :
   xMapping/src/xMapping.cc
   xMapping/src/xMapping.h

 and
   xFEM/src/xRegularGrid.cc
   modified to reflect xMapping::interiorCheck interface change
parent 6371e2c6
......@@ -146,9 +146,9 @@ int xRegularGrid ::GetElementsForPoint(const xtensor::xPoint &p, std::set<mEntit
assert(loc != nullptr);
if (e.find(loc) == e.end())
{
double u, v, w;
xPoint uvw;
xmapping::xMapping *mapping = xMappingBuilderHolderSingleton::instance().buildMapping(*loc);
if (mapping->interiorCheck(p, u, v, w))
if (mapping->interiorCheck(p, uvw))
{
e.insert(loc);
nbr++;
......
......@@ -236,6 +236,28 @@ vector3d compute_normal_to_solid_face(const xMapping &solidmapping, const unsign
xMapping::xMapping(xReferenceElementType _type) : type(_type) {}
xtensor::xPoint xMapping::eval(const xtensor::xPoint &uvw) const
{
xtensor::xPoint xyz;
eval(uvw[0], uvw[1], uvw[2], xyz[0], xyz[1], xyz[2]);
return xyz;
}
tensor xMapping::deval(const xtensor::xPoint &uvw) const
{
tensor dxdu;
deval(uvw[0], uvw[1], uvw[2], dxdu(0, 0), dxdu(1, 0), dxdu(2, 0), dxdu(0, 1), dxdu(1, 1), dxdu(2, 1), dxdu(0, 2), dxdu(1, 2),
dxdu(2, 2));
return dxdu;
}
xtensor::xBoundingBox xMapping::boundingBox() const
{
xtensor::xBoundingBox bb;
boundingBox(bb.min, bb.max);
return bb;
}
double xMapping::pushBack(double u, double v, double w, size_t vsize, vector3d *gr) const
{
tensor jInv;
......@@ -244,6 +266,11 @@ double xMapping::pushBack(double u, double v, double w, size_t vsize, vector3d *
return detJac;
}
double xMapping::pushBack(const xtensor::xPoint &uvw, size_t vsize, vector3d *vec) const
{
return pushBack(uvw[0], uvw[1], uvw[2], vsize, vec);
}
double xMapping::pushBack(double u, double v, double w, size_t vsize, tensor *gr) const
{
tensor jInv;
......@@ -265,6 +292,11 @@ double xMapping::pushBack(double u, double v, double w, size_t vsize, tensor *gr
return detJac;
}
double xMapping::pushBack(const xtensor::xPoint &uvw, size_t vsize, tensor *tens) const
{
return pushBack(uvw[0], uvw[1], uvw[2], vsize, tens);
}
bool xMapping::invert(double xp, double yp, double zp, double &Upos, double &Vpos, double &Wpos) const
{
#define NR_PRECISION 1.e-6
......@@ -298,6 +330,11 @@ bool xMapping::invert(double xp, double yp, double zp, double &Upos, double &Vpo
return true;
}
bool xMapping::invert(const xtensor::xPoint &xyz, xtensor::xPoint &uvw) const
{
return invert(xyz[0], xyz[1], xyz[2], uvw[0], uvw[1], uvw[2]);
}
// certainly not the fastest way !!
double xMapping::detJac(double u, double v, double w) const
{
......@@ -305,6 +342,12 @@ double xMapping::detJac(double u, double v, double w) const
return jacInverse(u, v, w, t);
}
double xMapping::detJac(const xtensor::xPoint &uvw) const
{
tensor t;
return jacInverse(uvw, t);
}
/**
This function is valid for any mesh mapping
*/
......@@ -391,7 +434,10 @@ double xMapping::jacInverse(double u, double v, double w, tensor &Invjac) const
Invjac(2, 2) = (jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) / DetJac;
return DetJac;
}
double xMapping::jacInverse(const xtensor::xPoint &uvw, tensor &Invjac) const
{
return jacInverse(uvw[0], uvw[1], uvw[2], Invjac);
}
bool xMapping::inReferenceElement(double u, double v, double w) const
{
const double eps = 1.e-6;
......@@ -427,18 +473,13 @@ bool xMapping::inReferenceElement(double u, double v, double w) const
return true;
}
bool xMapping::interiorCheck(const xtensor::xPoint &p, double &u, double &v, double &w) const
{
xtensor::xPoint pMin, pMax;
const double eps = 1.e-5;
boundingBox(pMin, pMax);
bool xMapping::inReferenceElement(const xtensor::xPoint &uvw) const { return inReferenceElement(uvw[0], uvw[1], uvw[2]); }
if (p(0) > pMax(0) + eps || p(1) > pMax(1) + eps || p(2) > pMax(2) + eps || p(0) < pMin(0) - eps || p(1) < pMin(1) - eps ||
p(2) < pMin(2) - eps)
return false;
if (!invert(p(0), p(1), p(2), u, v, w)) return false;
if (!inReferenceElement(u, v, w)) return false;
bool xMapping::interiorCheck(const xtensor::xPoint &p, xtensor::xPoint &uvw) const
{
if (!(boundingBox().contains(p))) return false;
if (!invert(p, uvw)) return false;
if (!inReferenceElement(uvw)) return false;
return true;
}
......@@ -461,13 +502,20 @@ void xMapping::COG(double &u, double &v, double &w) const
break;
default:
{
std::cout << "Eroor : COG unknown for type " << static_cast<int>(type) << " in " << __FILE__ << ":" << __LINE__
std::cout << "Error : COG unknown for type " << static_cast<int>(type) << " in " << __FILE__ << ":" << __LINE__
<< std::endl;
break; // need to throw an exception here
}
}
}
xtensor::xPoint xMapping::COG() const
{
xtensor::xPoint uvw;
COG(uvw[0], uvw[1], uvw[2]);
return uvw;
}
vector3d xMapping::normalVector(const unsigned short *face_vertices, double s, double t) const
{
switch (type)
......
......@@ -10,6 +10,7 @@
#include <vector>
// xtensor
#include "xBoundingBox.h"
#include "xPoint.h"
#include "xTensor2.h"
#include "xVector.h"
......@@ -21,62 +22,103 @@ namespace xmapping
using tensor = xtensor::xTensor2<double>;
using vector3d = xtensor::xVector<double>;
// \note xReferenceElementType is defined by AOMD_Internals. we need to keep it this way since Integrator from AOMD uses it.
// We will could have our own when we have our own xIntegrator Class.
/// xMapping define a mapping from a reference domain to the space domain
/// xMapping is virtual class to define a mapping from a reference domain to the space domain
/*!
* \brief The xMapping class define a mapping from a reference domain to the space domain.
* u, v, w refere to coordinate in the reference space, x, y, z refer to coordinate in the physical space.
* derived class must implement a minima the following pur virtual function :
* virtual void eval(double u, double v, double w, double &x, double &y, double &z) const = 0;
* virtual void deval(double u, double v, double w, double &dxdu, double &dydu, double &dzdu, double &dxdv, double &dydv,
* double &dzdv, double &dxdw, double &dydw, double &dzdw) const = 0;
* virtual void boundingBox(xtensor::xPoint &min, xtensor::xPoint &max) const = 0;
* virtual int order() const = 0;
* virtual int geomOrder() const = 0;
*/
class xMapping
{
public:
xMapping(xReferenceElementType type);
virtual ~xMapping() = default;
/// return the type of the refence element
inline xReferenceElementType getType() const { return type; }
//! Inversion of the mapping by a newton algorithm
/// return the coordinates x,y,z in real space corresponding to the local cordinate u,v,w
virtual void eval(double u, double v, double w, double &x, double &y, double &z) const = 0;
/// return the coordinates xyz in real space corresponding to the local cordinate uvw
virtual xtensor::xPoint eval(const xtensor::xPoint &uvw) const;
/// return the partial derivatives of real coordinates xyz with regard to uvw at local cordinate u,v,w
virtual void deval(double u, double v, double w, double &dxdu, double &dydu, double &dzdu, double &dxdv, double &dydv,
double &dzdv, double &dxdw, double &dydw, double &dzdw) const = 0;
/// return the partial derivatives of real coordinates xyz with regard to uvw at local cordinate u,v,w in a tensor such as
//! (i,j) = dx_i/du_j
virtual tensor deval(const xtensor::xPoint &uvw) const;
/// Inversion of the mapping by a newton algorithm. return true if succcessful, false other wise
virtual bool invert(double x, double y, double z, double &u, double &v, double &w) const;
/// Inversion of the mapping by a newton algorithm. return true if succcessful, false other wise
virtual bool invert(const xtensor::xPoint &xyz, xtensor::xPoint &uvw) const;
/*!
computation of the inverse of the jacobian at u,v,w (invjac),
where jacobian = d(x,y,z)/d(u,v,w) is a 3x3 matrix
returns the determinant of the jacobian (not of the inverse).
*/
virtual double jacInverse(double u, double v, double w, tensor &invjac) const;
/*!
Jacobian = d(x,y,z)/d(u,v,w) is a 3x3 matrix
computation of the inverse of the jacobian at (u,v,w),
computation of the inverse of the jacobian at point uvw (invjac),
where jacobian = d(x,y,z)/d(u,v,w) is a 3x3 matrix
returns the determinant of the jacobian (not of the inverse).
*/
virtual double jacInverse(double x, double y, double z, tensor &) const;
virtual double jacInverse(const xtensor::xPoint &uvw, tensor &invjac) const;
//! returns determinant of jacobian at (u,v,w)
virtual double detJac(double u, double v, double w) const;
//! returns determinant of jacobian at point uvw
virtual double detJac(const xtensor::xPoint &uvw) const;
/**
tells if a point (u,v,w) is inside the reference element or not
For EDGE, return true if ( -1.-eps <= u <= 1.+eps(
for TRI, return true if (u >= -eps && v >= -eps && 1.-u-v => -eps)
for TRI, return true if (u >= - && v >= -eps && 1.-u-v => -eps)
for QUAD, return true if (u >= (-1.-eps) && u <= (1+eps) && v >= (-1.0-eps) && v <= (1.0+eps))
for TET, return true if (u => (-eps && v => -eps && w => -eps && 1.-u-v-w => -eps)
for HEX, return true if ( (1+eps) => u, v, w => (-1.-eps)==)
for PRISM, return true if (-1-eps <= w <= 1+eps && u,v >= -eps && 1.-u-v >= -eps )
*/
virtual bool inReferenceElement(double u, double v, double w) const;
virtual bool inReferenceElement(const xtensor::xPoint &uvw) const;
/**
checks if a point in REAL coordinates is inside the element ent
if it's the case, returns local coordinates
*/
virtual bool interiorCheck(const xtensor::xPoint &p, double &u, double &v, double &w) const;
/// returns the center of gravity of the reference element in local coordinates
/*!
* (it is not in general the local coordinate of the center of gravity of the element ... even if it is the case for some
* linear mapping)
*/
virtual bool interiorCheck(const xtensor::xPoint &p, xtensor::xPoint &uvw) const;
/// returns the local coordinates of the COG of the reference element
virtual void COG(double &u, double &v, double &w) const;
/// Computes gradients in real world knowing them in reference world, vector version
/// returns the local coordinates of the COG of the reference element
virtual xtensor::xPoint COG() const;
/// Computes gradients at local cordinate u,v,w real world knowing them in reference world, vector version
//! vsize is the number of vector to push,
//! on entry vec is a pointer to contiguous vector3d, supposed to represent partialderivative with regrad to u,v,w.
//! on exit vec is a pointer to contiguous vector3d, supposed to represent gradient in x,y,z.
//! the returned value is determinent of the jacobian
//! \note could use span when switching to c++20
virtual double pushBack(double u, double v, double w, size_t vsize, vector3d *vec) const;
/// Computes gradients in real world knowing them in reference world, tensor version
virtual double pushBack(double u, double v, double w, size_t vsize, tensor *tens) const;
/// Computes gradients at local cordinate point uvw real world knowing them in reference world, vector version
//! vsize is the number of vector to push,
//! on entry vec is a pointer to contiguous vector3d, supposed to represent partialderivative with regrad to u,v,w.
//! on exit vec is a pointer to contiguous vector3d, supposed to represent gradient in x,y,z.
//! the returned value is determinent of the jacobian
virtual double pushBack(const xtensor::xPoint &uvw, size_t vsize, vector3d *vec) const;
inline double pushBack(double u, double v, double w, std::vector<vector3d> &vec) const
{
return pushBack(u, v, w, vec.size(), &(*vec.begin()));
return pushBack(u, v, w, vec.size(), vec.data());
}
inline double pushBack(const xtensor::xPoint &uvw, std::vector<vector3d> &vec) const
{
return pushBack(uvw, vec.size(), vec.data());
}
/// Computes gradients in real world knowing them in reference world, tensor version
virtual double pushBack(double u, double v, double w, size_t vsize, tensor *tens) const;
/// Computes gradients in real world knowing them in reference world, tensor version
virtual double pushBack(const xtensor::xPoint &uvw, size_t vsize, tensor *tens) const;
/// compute the exterior normal to the face face_vertices at point s,t of the face
/*!
......@@ -106,28 +148,14 @@ class xMapping
* it also work with line, making the assumption that the normal is orthogonal to the third axis (dxdv set to e3).
*/
virtual vector3d normalVector(double u, double v, double w) const;
//---------------------------------------------//
// P u r e V i r t u a l M e m b e r s //
//---------------------------------------------//
/**
eval and deval functions have to be re-written for different mesh
mappings (curvilinear, bezier,...).
eval computes x(u,v,w) , y(u,v,w) and z(u,v,w)
deval computes derivatives dx/du , dy/du ...
*/
virtual void eval(double u, double v, double w, double &x, double &y, double &z) const = 0;
virtual void deval(double u, double v, double w, double &dxdu, double &dydu, double &dzdu, double &dxdv, double &dydv,
double &dzdv, double &dxdw, double &dydw, double &dzdw) const = 0;
/**
computes the bounding box of a mesh entity
*/
/// computes the bounding box of the mapping in real space.
virtual void boundingBox(xtensor::xPoint &min, xtensor::xPoint &max) const = 0;
/// computes the bounding box of the mapping in real space.
virtual xtensor::xBoundingBox boundingBox() const;
/// return the polynomial order of the mapping (if it make sense ...), used for chosing integration rule
virtual int order() const = 0;
/// return the polynomial order of the mapping (if it make sense ...), used for chosing integration rule
//! similar to order. probably here for historical reason ... \note I need to figure this out
virtual int geomOrder() const = 0;
protected:
......
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