Commit b71ddff8 authored by Nicolas CHEVAUGEON's avatar Nicolas CHEVAUGEON

[xFEM, xMapping] xGeomElem code is simplifyed mostlty by using improved xMapping interface.

xGeomElem code is simplifyed (and improved ?) in the following way :

- All the members now have version with the correct camelCase style. The old version are kept (for example GetDetJac is "cloned" to getDetJac), but are marked as "deprecated" in the their doxygen documentation. I tryed to use the [[deprecated]] attribute on them but unfortunatly, it's c++14 and we still can't use this version of the standard. I kept all the old version there since xGeomElem is used all over the places, but I would like to remove them in a near future ...

- xBoundingBox xGeomElem::getBoundingBox() const take the place of void xGeomElem::GetBoundingBox(xtensor::xPoint& min, xtensor::xPoint& max) const.

- The implementation take advantages of improved interface to xMapping.h
  for example :
  xtensor::xPoint xGeomElem ::getCDGxyz() const
  {
   double u, v, w, x, y, z;
   mapping->COG(u, v, w);
   mapping->eval(u, v, w, x, y, z);
   return xtensor::xPoint(x, y, z);
  }
  is replaced by
  xtensor::xPoint xGeomElem::getCDGxyz() const { return mapping->eval(mapping->COG()); }
  and
  void xGeomElem::setUVWForVertex(int inod) which was almost 200 lines
  is replaced by
  void xGeomElem ::setUVWForVertex(int inod) { uvw = mapping->getMappingVertex(inod); }
  where mapping->getMappingVertex(inod); just return data that where already in xMapping.

- setting/clearing  the mapping or integrator members of the class was complicated because we have to keep track of differents life  time for  these pointers depending if they where defaulty allocated or pass by the user and eventulaly shared by differents
  xGeomElem. This is now simplifyed using a unique_ptr and classic ptr, allowing to use the implicitly declared destructor for
  xGEomElem.
  This pave the way for better mapping management in the  futur (I'm thinking of high order mapping that would be expensive
  to construct, that we might wan't to store in the builder ...)

- we now have move constructor, move assignment and swap for xGeomElem, making them more usable in standard containers.

Note on xGeomElem :
 - I think set/getOrientation are not very usefull, probably bug prone, and I think they should be removed.
 - I'm not sure we still need the integrator to be included in the xGeomElem, Integration behing performed in most cases using an
   xIntegrationRule  and xCommandOnGeomElem that could set directly the current point instead of setting it using both xGeomElem::setIntegrationPointNumberForDegree and xGeomElem::setUVW(int ipt)

4 access functions have been added in class xMapping to access data that where already there, but not accessible outside xMapping.cc :
   -size_t getNbMappingVertices() const;
   -xtensor::xPoint getMappingVertex(size_t i) const;
   -static size_t getNbMappingVertices(xReferenceElementType type);
   -static xtensor::xPoint getMappingVertex(xReferenceElementType type, size_t i);
parent a7cd5a68
This diff is collapsed.
/*
/*
This file is a part of eXlibris C++ Library
under the GNU Lesser General Public License.
See the NOTICE.md & LICENSE.md files for terms
......@@ -7,6 +7,7 @@
#ifndef _xGEOMETRICAL_ELEMENT_H
#define _xGEOMETRICAL_ELEMENT_H
#include <memory>
// aomd trellis
#include "mEntity.h"
// xquadrature
......@@ -23,76 +24,131 @@ class xMapping;
namespace xfem
{
/// xGeomElem combines an entity, a xMapping, an xIntegrator and a current point it the reference configuration.
///! It is heavily used to evaluate functions inside elements either in reference space or in geometrical space at current point.
///! or push gradient of function with regard to localcoordinate to global coordinate via the mapping.
///! UVW refers to the current point in the reference configuration. By default uvw is set to {0.,0.,0.}
///! UVW can be set manually with setUVW(const xPoint &) or setUVWforXYZ(const xPoint &)
///! or it can be set using the xIntegrator by setting the ith point in the current integration rule setUVW(int)
class xGeomElem
{
public:
xGeomElem(const xGeomElem& other);
xGeomElem& operator=(const xGeomElem& other);
xGeomElem(AOMD::mEntity*); // Test K ¤ 12/02/09
/// Construct a GeomElem from a mEntity.
//! The mapping is constructed using xMappingBuilderHolderSingleton::instance().buildMapping(*e)
//! And the integrator is constructed by calling xquadrature::xGaussIntegrator.
//! The mapping and the integrator are detroyed with the GeomElem or when calling setMapping
//! and setIntegrator respectively
xGeomElem(AOMD::mEntity* e);
/// Construct a GeomElem, from a mapping, and an integrator
//! xGeomElem don't take ownership of the mapping nor the integrator in this case.
//! They must not be destroyed before they are reset or the xGeomElem is destroy.
xGeomElem(AOMD::mEntity* e, xmapping::xMapping* m, xquadrature::xIntegrator* i);
~xGeomElem();
//
/// Copy constructor. Note that the mapping and the integrator are fully copied if they where owned by other.
xGeomElem(const xGeomElem& other);
xGeomElem(xGeomElem&& other) noexcept;
xGeomElem& operator=(xGeomElem other);
friend void swap(xGeomElem& l, xGeomElem& r) noexcept;
/// return the pointer to the mesh entity.
AOMD::mEntity* getEntity() const { return const_cast<AOMD::mEntity*>(pent); }
xtensor::xPoint getXYZ() const;
/// return the local coordinates of the current point.
const xtensor::xPoint& getUVW() const { return uvw; }
xtensor::xPoint getCDGxyz() const;
/// return the global coordinates of the curreent point (using the mapping)
xtensor::xPoint getXYZ() const;
/// return the center of gravity of the ref element in local coordinate.
xtensor::xPoint getCDGuvw() const;
// computes volume (3D), area (2D) or length (1D).
/// return the center of gravity of the ref element pushed to global coordinate.
xtensor::xPoint getCDGxyz() const;
/// computes volume (3D), area (2D) or length (1D), using the current integrator, with rule of order 0
double getMeasure();
/// return the exterior normal to the GeomElem, on the given boder.
//! The border must one of the dim-1 element that bound the entitie's geom Elem or throw.
xtensor::xVector<> normalVector(const AOMD::mEntity& pborder) const;
/// [deprecated("Use xtensor::xVector<> normalVector(const mEntity& pborder) const"]]
void normalVector(AOMD::mEntity* border, xtensor::xVector<>& n) const;
/// return the bounding box of the mapping.
xtensor::xBoundingBox getBoundingBox() const;
/// [[deprecated("Use getBoundingBox() instead")]]
void GetBoundingBox(xtensor::xPoint& min, xtensor::xPoint& max) const;
bool PointInElement(xtensor::xPoint& xyz) const;
/// return true if point xyz is inside the element (according to the mapping)
bool pointInElement(const xtensor::xPoint& xyz) const;
// [[deprecated("Use bool PointInElement(const xtensor::xPoint& xyz) const; instead")]]
bool PointInElement(const xtensor::xPoint& xyz) const;
double getDetJac() const;
/// [[deprecated("Use getDetJac() instead"]]
double GetDetJac() const;
const double& GetWeight() const { return Weight; }
unsigned int GetNbIntegrationPoints() const { return NbIntegrationPoints; }
/// return the Weight of the current point, as set by the integration rule.
const double& getWeight() const { return Weight; }
/// [[deprecated("Use getWeight() instead"]]
const double& GetWeight() const { return getWeight(); }
unsigned int getNbIntegrationPoints() const { return NbIntegrationPoints; }
/// [[deprecated("Use getNbIntegrationPoints() t instead"]]
unsigned int GetNbIntegrationPoints() const { return getNbIntegrationPoints(); }
/// get/set orientation of the GeomElem : 1 : same as the entity, -1 : reversed. the jacobian is multiply by the orientation
/// when computing detJac. \note theses function seems quite useless and impose some book keeping. We probably should remove
/// them.
int getOrientation() const { return Ori; }
void setOrientation(int s) { Ori = (s >= 0) ? 1 : -1; }
/// [[deprecated("Use getOrientation instead"]]
int GetOrientation() const { return Ori; }
void SetOrientation(int s)
{
if (s >= 0.)
Ori = 1;
else
Ori = -1;
}
// to select the integration rule
/// [[deprecated("Use setOrientation instead"]]
void SetOrientation(int s) { Ori = (s >= 0) ? 1 : -1; }
/// set the degree of the integrator.
void setIntegrationPointNumberForDegree(const int& degtot);
/// [[deprecated("Use setIntegrationPointNumberForDegree(const int& degtot) instead")]]
void SetIntegrationPointNumberForDegree(const int& degtot);
unsigned int getOrder() { return order; } // Give the order (or degree)
/// return the order of the integration rule (as set by ) setIntegrationPointNumberForDegree
unsigned int getOrder() { return order; }
/// set the current local point to the ipt point of the integration rule.
//! ipt must be lower than the number of points in the current integration rule (getNbIntegrationPoints());
void setUVW(int ipt);
/// set the current point local coordinates
void setUVW(const xtensor::xPoint& uvw_);
/// set the current point by inverting the mapping at xyz. Can fail if xyz is not on the element.
void setUVWForXYZ(const xtensor::xPoint& xyz);
/// set UVW at the mapping vertex inod (see xMapping)
void setUVWForVertex(int inod);
/// return the current point number in the integration rule.
int getCurrentIntegrationPointId() const { return CurrentIntegrationPoint; }
xtensor::xVector<>& pushBack(xtensor::xVector<>& inx) const;
xtensor::xTensor2<>& pushBackRight(xtensor::xTensor2<>& inx) const;
xtensor::xTensor2<>& pushBackRightTranspose(xtensor::xTensor2<>& inx) const;
/// [[deprecated("Use pushBack(xtensor::xVector<>& inx) instead")]]
xtensor::xVector<>& PushBack(xtensor::xVector<>& inx) const;
/// [[deprecated("Use pushBackRight(xtensor::xVector<>& inx) instead")]]
xtensor::xTensor2<>& PushBackRight(xtensor::xTensor2<>& inx) const;
/// [[deprecated("Use pushBackRightTranspose(xtensor::xVector<>& inx) instead")]]
xtensor::xTensor2<>& PushBackRightTranspose(xtensor::xTensor2<>& inx) const;
/// return the element's zone tag.
int getZone() const;
/// return the mapping
const xmapping::xMapping* getMapping() const { return mapping; }
/// set the mapping. if mapping set this way, the user is responsible for the pointed to lifetime.
void setMapping(xmapping::xMapping* m);
/// return the quadrature.
const xquadrature::xIntegrator* getIntegrator() const { return integrator; }
/// set the integrator. if setted this way, the user is responsible for the life time.
void setIntegrator(xquadrature::xIntegrator* integrator);
private:
// \note : need to be a pointer to cope with operator =.
const AOMD::mEntity* pent;
int Ori;
bool default_constructor;
bool default_Integrator;
bool default_Mapping;
xmapping::xMapping* mapping;
xquadrature::xIntegrator* integration_points;
int order;
const AOMD::mEntity* pent = nullptr;
int Ori = 1;
// \note the mapping and the integrator can be owned or not by the xGeomElem.
// so we have a unique pointer ownedmapping set to the mapping if xGeomElem owns it,
// and mapping point to the ownedmapping target in this case.
// if xGeomElem don't own it, ownedmapping is set to nullptr and mapping point to what ever the user want,
// but wat is pointed to won't be destroyed at detruction of xGeomElem.
// the same idiom is used for xintegrator.
std::unique_ptr<xmapping::xMapping> ownedmapping = nullptr;
xmapping::xMapping* mapping = nullptr;
std::unique_ptr<xquadrature::xIntegrator> ownedintegrator = nullptr;
xquadrature::xIntegrator* integrator = nullptr;
int order = 0;
xtensor::xPoint uvw;
int CurrentIntegrationPoint;
unsigned int NbIntegrationPoints;
double Weight;
int CurrentIntegrationPoint = 0;
unsigned int NbIntegrationPoints = 1;
double Weight = 1.;
xGeomElem() = default;
};
} // namespace xfem
......
......@@ -12,27 +12,26 @@
namespace xmapping
{
const double edge_vert_pos[2] = {-1., 1.};
const vector3d vertex_normal_of_edge[2] = {{-1., 0., 0.}, {1., 0., 0.}};
const double tri_vert_pos[3][2] = {{0., 0.}, {1., 0.}, {0., 1.}};
const double quad_vert_pos[4][2] = {{-1., -1.}, {1., -1.}, {1., 1.}, {-1., 1.}};
const double tet_vert_pos[4][3] = {{0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}};
const double prism_vert_pos[8][3] = {{0., 0., -1.}, {1., 0., -1.}, {0., 1., -1.}, {0., 0., 1.}, {1., 0., 1.}, {0., 1., 1.}};
const double hex_vert_pos[8][3] = {{-1., -1., -1.}, {1., -1., -1.}, {1., 1., -1.}, {-1., 1., -1.},
{-1., -1., 1.}, {1., -1., 1.}, {1., 1., 1.}, {-1., 1., 1.}};
const vector3d vertex_normal_of_edge[2] = {{-1., 0., 0.}, {1., 0., 0.}};
const unsigned short edge_vertex_of_tri[3][2] = {{0, 1}, {1, 2}, {2, 0}};
const vector3d edge_normal_of_tri[3] = {{0., -1., 0.}, {1. / sqrt(2.), 1. / sqrt(2.), 0.}, {-1., 0., 0.}};
const unsigned short swapcase_tri[6][3] = {{0, 1, 2}, {2, 0, 1}, {1, 2, 0}, {0, 2, 1}, {2, 1, 0}, {1, 0, 2}};
const double quad_vert_pos[4][2] = {{-1., -1.}, {1., -1.}, {1., 1.}, {-1., 1.}};
const unsigned short edge_vertex_of_quad[4][2] = {{0, 1}, {1, 2}, {3, 2}, {0, 3}};
const vector3d edge_normal_of_quad[4] = {{0., -1., 0.}, {1., 0., 0.}, {0., 1., 0.}, {-1., 0., 0.}};
const unsigned short swapcase_quad[8][4] = {{0, 1, 2, 3}, {3, 0, 1, 2}, {2, 3, 0, 1}, {1, 2, 3, 0},
{1, 0, 3, 2}, {2, 1, 0, 3}, {3, 2, 1, 0}, {0, 3, 2, 1}};
const double tet_vert_pos[4][3] = {{0., 0., 0.}, {1., 0., 0.}, {0., 1., 0.}, {0., 0., 1.}};
const unsigned short face_vertex_of_tet[4][3] = {{0, 1, 2}, {0, 1, 3}, {1, 2, 3}, {0, 2, 3}};
const vector3d face_normal_of_tet[4] = {
{0., 0., -1.}, {0., -1., 0.}, {1. / sqrt(3.), 1. / sqrt(3.), 1. / sqrt(3.)}, {-1., 0., 0.}};
const double hex_vert_pos[8][3] = {
{-1., -1., -1.}, {1., -1., -1.}, {1., 1., -1.}, {-1., 1., -1.}, {-1., -1., 1.}, {1., -1., 1.}, {1., 1., 1.}, {-1., 1., 1.},
};
const unsigned short face_vertex_of_hex[6][4] = {{0, 1, 2, 3}, {0, 1, 5, 4}, {1, 2, 6, 5},
{3, 2, 6, 7}, {0, 3, 7, 4}, {4, 5, 6, 7}};
const vector3d face_normal_of_hex[6] = {{0., 0., -1.}, {0., -1., 0.}, {1., 0., 0.}, {0., 1., 0.}, {-1., 0., 0.}, {0., 0., 1.}};
......@@ -605,4 +604,62 @@ vector3d xMapping::normalVector(double u, double v, double w) const
}
}
size_t xMapping::getNbMappingVertices(xReferenceElementType type)
{
switch (type)
{
case xReferenceElementType::VERTEX:
return 1;
case xReferenceElementType::EDGE:
return 2;
case xReferenceElementType::TRI:
return 3;
case xReferenceElementType::QUAD:
return 4;
case xReferenceElementType::TET:
return 4;
case xReferenceElementType::PYRAMID:
return 5;
case xReferenceElementType::PRISM:
return 6;
case xReferenceElementType::HEX:
return 8;
default:
throw;
}
}
xtensor::xPoint xMapping::getMappingVertex(xReferenceElementType type, size_t i)
{
switch (type)
{
case xReferenceElementType::VERTEX:
assert(i < 1);
return xtensor::xPoint{0., 0., 0.};
case xReferenceElementType::EDGE:
assert(i < 2);
return xtensor::xPoint{edge_vert_pos[i], 0., 0.};
case xReferenceElementType::TRI:
assert(i < 3);
return xtensor::xPoint{tri_vert_pos[i][0], tri_vert_pos[i][1], 0.};
case xReferenceElementType::QUAD:
assert(i < 4);
return xtensor::xPoint{quad_vert_pos[i][0], quad_vert_pos[i][1], 0.};
case xReferenceElementType::TET:
assert(i < 4);
return xtensor::xPoint{tet_vert_pos[i][0], tet_vert_pos[i][1], tet_vert_pos[i][2]};
case xReferenceElementType::PRISM:
assert(i < 6);
return xtensor::xPoint{prism_vert_pos[i][0], prism_vert_pos[i][1], prism_vert_pos[i][2]};
case xReferenceElementType::HEX:
assert(i < 8);
return xtensor::xPoint{hex_vert_pos[i][0], hex_vert_pos[i][1], hex_vert_pos[i][2]};
default:
throw;
}
}
xtensor::xPoint xMapping::getMappingVertex(size_t i) const { return getMappingVertex(type, i); }
size_t xMapping::getNbMappingVertices() const { return getNbMappingVertices(type); }
} // namespace xmapping
......@@ -157,6 +157,10 @@ class xMapping
/// 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;
size_t getNbMappingVertices() const;
xtensor::xPoint getMappingVertex(size_t i) const;
static size_t getNbMappingVertices(xReferenceElementType type);
static xtensor::xPoint getMappingVertex(xReferenceElementType type, size_t i);
protected:
xReferenceElementType type;
......
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