Commit 22fa74d5 authored by Nicolas CHEVAUGEON's avatar Nicolas CHEVAUGEON

[XFEM, xExport, xInterface/xOctree] dependences to AOMD's octree is removed and related changes.

xFEM/src/xOctreeGrid.h is a new file. It introduce a octree grid search
  structure to use in place of the AOMD::octree that was used in xMesh
  for the locateOctree members functions. It permits to construct an octree
  search structure on top of a collection of mesh entity.

xFEM/src/xRegularGrid.h/cc has been modifyed so that :
 - it uses the xBoundingBox class.
 - The member function name have been changed to reflect the camel case
   convention, common to the rest of the xFEM library.
 - getElementsForPoint now return a vector of {element, localcoordinate} pair,
   so that the local coordinate of the target point in each element that
   contains it is now accessible.
 - The code have been clarifyed and simplifyed.

xBoundingBox.h :
  - added 2 member function in xBoundingBox :
   /// return the center of the bounding box
   xPoint center() const { return (min + max) * 0.5; }
   /// return a vector from min to max
   xVector<double> diag() const { return xVector<double>{min, max}; }
  - Note that following a discussion with A. Salzmann, xBoundingBox could
    be moved in a near futur to namespace xgeom and the file might be part
    of a new library (xGeomTools maybe)

xFEM/src/xMesh.h/cc has been modifyed so that :
  - It uses the xBoundingBox class were usefull, implying some interface
    change :
     * std::pair<xPoint, xPoint> compute_bounding_box(mEntity* e)
       is changed to :
       xBoundingBox compute_bounding_box(const mEntity& e);
     * void compute_bounding_box(xPoint& min,xPoint& max) const;
       is "deprecated", but is still there during the transition.
       it is replaced by :
       xtensor::xBoundingBox compute_bounding_box() const;
  - locateElementOctree and locateElement have changed interface,
    reflecting the changes in xRegularGrid and the use of xOctreeGrid instead of
    AOMD's octree.
  - Better control of the search structure in xMesh, since now the octree or
    the regular grid can both be cleared, updated or accessed by calling
    respectively xMesh Members   createOctree(), clearOctree(), getOctreeGrid(),
    createGrid(), clearGrid(), getRegularGrid.
  - The search structure pointers are now managed in a unique_ptr.
  - xMesh::copyMeshInternal now copy the search structures if needed instead
    of throwing.
  - AOMD_SharedInfo.h, autopack.h, OctreeCreate2.h (from AOMD) are not
    included anymore.

xCrack/src/CrackPostpro.cc,  xCrack/src/ExportSensors.cc/h, modified so that it
  reflect the change in xMesh::locateElement

 xExport/src/xSensors.cc/h :
  -modified so that it reflect the change in xMesh::locateElement,
  -Trellis_Util::mPoint replaced by xtensor::xPoint
  -small changes to simplify the implementation
  -clang-format first the first time I think.

xInterface/xOctree/src/surf2LevelSet.cc/h
  -modifyed to use xBoundingBox
  -surface2LevelSet::getBB(xPoint&, xPoint&) is replaced by
   const xBoundingBox &getBB() const { return BB; }
  -surface2LevelSet members xPoint BBmin, BBmax are replaced xBoundingBox BB;
  -a bug was found in the process :
   In surface2LevelSet::createMeshBBox that create an xMesh out of grid,
   the grid point coordinate were computed as :
    double xyz[3] = {BB.min(0) + i * elemSizeX, BB.min(1) +
                                 j * elemSizeY, BB.min(2) + k * elemSizeY};
   instead of :
     double xyz[3] = {BB.min(0) + i * elemSizeX, BB.min(1) +
                                  j * elemSizeY, BB.min(2) + k * elemSizeZ};
   The bug is corrected, but reference of test case
   xinterface-xoctree-test/surface2LevelSet had to be changed.

All modifyed files have been clang-fomated.
parent 661796b4
......@@ -34,7 +34,6 @@
using namespace std;
using namespace xfem;
using namespace std;
using namespace AOMD;
namespace xcrack
......@@ -284,20 +283,15 @@ void lCrack::extractSIFs(const xField<>& disp_l, std::string filename, bool auto
JEM.setUVW(i); // Upos, Vpos, Wpos and JacMatrix are set
xtensor::xPoint xyz = JEM.getXYZ();
if (debug) printf("XPos YPos ZPos %12.5e %12.5e %12.5e\n", xyz(0), xyz(1), xyz(2));
std::set<mEntity*> eset;
mesh->locateElement(xyz, eset); // todo
if (debug) cout << "eset size is " << eset.size() << endl;
if (eset.size() == 0)
{
// cout << "warning xyz out of the box " << xyz << endl;
out_of_the_box = true;
// assert(0);
}
double eset_size = (double)eset.size();
for (std::set<mEntity*>::const_iterator its = eset.begin(); its != eset.end(); ++its)
auto list_pe_uvw = mesh->locateElement(xyz); // todo
if (list_pe_uvw.empty()) out_of_the_box = true;
double eset_size = (double)std::distance(list_pe_uvw.begin(), list_pe_uvw.end());
for (auto pe_uvw : list_pe_uvw)
{
mEntity* e = *its; // IxElement of the FE mesh
xfem::xGeomElem geom_appro(e);
const AOMD::mEntity* pe = pe_uvw.first;
xfem::xGeomElem geom_appro(const_cast<AOMD::mEntity*>(pe));
xMaterial* mat = xMaterialManagerSingleton::instance().getMaterial(&geom_appro);
const xTensors* properties = mat->getProperties();
MatInfo[0] = properties->scalar("YOUNG_MODULUS");
......@@ -313,7 +307,7 @@ void lCrack::extractSIFs(const xField<>& disp_l, std::string filename, bool auto
dalpha = gradqloc * axis;
double alpha = GetqLoc(postid, xyz, JBoxDimension, JFrontPoint, JMap);
if (debug) cout << "alpha is " << alpha << endl;
geom_appro.setUVWForXYZ(JEM.getXYZ());
geom_appro.setUVW(pe_uvw.second);
double Jhh_l;
xtensor::xVector<> Ih_l;
......
......@@ -23,27 +23,24 @@ void xeExportSensors::setMesh(xMesh* mesh)
}
sensor_to_elt.clear();
std::map<int, xtensor::xPoint>::const_iterator it = sensors.begin();
int i = 0;
for (; it != sensors.end(); ++it, ++i)
for (const auto& id_pt : sensors)
{
// int is = it->first;
xtensor::xPoint p = it->second;
std::set<mEntity*> elts;
mesh->locateElement(p, elts);
if (elts.size() == 0)
xtensor::xPoint p = id_pt.second;
auto elt_uvw_list = mesh->locateElement(p);
if (elt_uvw_list.empty())
{
cerr << "the sensor located at " << p << " is not in the mesh " << endl;
abort();
}
std::set<mEntity*>::const_iterator it_elts = elts.begin();
bool check = false;
for (; it_elts != elts.end(); ++it_elts)
for (auto elt_uvw : elt_uvw_list)
{
if (test(*it_elts))
AOMD::mEntity* pe = const_cast<AOMD::mEntity*>(elt_uvw.first);
if (test(pe))
{
xfem::xGeomElem* geo = new xfem::xGeomElem(*it_elts);
geo->setUVWForXYZ(p);
xfem::xGeomElem* geo = new xfem::xGeomElem(pe);
geo->setUVW(elt_uvw.second);
if (debug) cout << " in SetMesh, location of sensor is " << geo->getXYZ() << endl;
sensor_to_elt.insert(make_pair(i, geo));
check = true;
......@@ -54,6 +51,7 @@ void xeExportSensors::setMesh(xMesh* mesh)
cerr << "the sensor located at " << p << " is is in the mesh but not in a valid element" << endl;
abort();
}
++i;
}
}
......
......@@ -38,11 +38,7 @@ class xeExportSensors
}
~xeExportSensors()
{
std::map<int, xfem::xGeomElem*>::iterator itstart = sensor_to_elt.begin();
for (; itstart != sensor_to_elt.end(); itstart++)
{
delete itstart->second;
}
for (auto& id_geom : sensor_to_elt) delete id_geom.second;
}
void setMesh(xMesh* m);
......
/*
/*
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
See the NOTICE.md & LICENSE.md files for terms
and conditions.
*/
#include <list>
#include <set>
// xfem
#include "xGeomElem.h"
#include "xMesh.h"
// xexport
#include "xSensors.h"
// AOMD
#include "mEntity.h"
// xtool
#include "workInProgress.h"
#include <list>
#include <set>
namespace xexport
{
void xSensors::registerSensor(const int nb_, const Trellis_Util::mPoint pt_, xfem::xMesh* mesh_)
{
point_container.insert(std::pair<int, Trellis_Util::mPoint>(nb_, pt_));
std::set<AOMD::mEntity*> elts;
mesh_->locateElement(pt_, elts);
// in // if all proc do have a elts.size()==0 then sensor is out of mesh otherwise it is just in some proc that are not
// this one. => reduce on max of elts.size to test
// if (max==0) => message
// else if (elts.size()==0) => nothing
// else treatement
assert(xtool::workInProgress());
if (elts.size()==0)
{
std::cout << "sensor is out of mesh, it is not set" << std::endl;
}
else
{
AOMD::mEntity* elt = *(elts.begin());
xfem::xGeomElem* geo = new xfem::xGeomElem(elt);
geo->setUVWForXYZ(pt_);
geo->SetIntegrationPointNumberForDegree(0);
geom_elem_container.insert(std::pair<int, xfem::xGeomElem*>(nb_, geo));
}
}
void xSensors::registerSensor(const int nb_, const xtensor::xPoint pt_, xfem::xMesh* mesh_)
{
point_container.insert(std::make_pair(nb_, pt_));
auto list_pe_uvw = mesh_->locateElement(pt_);
// in // if all proc do have a elts.size()==0 then sensor is out of mesh otherwise it is just in some proc that are not
// this one. => reduce on max of elts.size to test
// if (max==0) => message
// else if (elts.size()==0) => nothing
// else treatement
assert(xtool::workInProgress());
if (list_pe_uvw.empty())
std::cout << "sensor is out of mesh, it is not set" << std::endl;
else
{
auto pe_uvw = list_pe_uvw.front();
std::unique_ptr<xfem::xGeomElem> geo(new xfem::xGeomElem(const_cast<AOMD::mEntity*>(pe_uvw.first)));
geo->setUVW(pe_uvw.second);
geo->SetIntegrationPointNumberForDegree(0);
geom_elem_container.insert(std::make_pair(nb_, std::move(geo)));
}
}
void xSensors::reinit(xfem::xMesh* mesh_)
{
clean();
std::map<int, Trellis_Util::mPoint>::const_iterator it=point_container.begin();
std::vector<int> to_remove;
for (; it!=point_container.end(); ++it)
void xSensors::reinit(xfem::xMesh* mesh_)
{
clean();
std::vector<int> to_remove;
for (const auto& id_point : point_container)
{
int nb = id_point.first;
const xtensor::xPoint& p = id_point.second;
auto list_pe_uvw = mesh_->locateElement(p);
if (list_pe_uvw.empty())
{
int nb = it->first;
Trellis_Util::mPoint p = it->second;
std::set<AOMD::mEntity*> elts;
mesh_->locateElement(p, elts);
if (elts.size()==0)
{
std::cout << "sensor " << nb << " is out of mesh, it is removed" << std::endl;
to_remove.push_back(nb);
}
else
{
AOMD::mEntity* elt = *(elts.begin());
xfem::xGeomElem* geo = new xfem::xGeomElem(elt);
geo->setUVWForXYZ(p);
geo->SetIntegrationPointNumberForDegree(0);
geom_elem_container.insert(std::pair<int, xfem::xGeomElem*>(nb,geo));
}
std::cout << "sensor " << nb << " is out of mesh, it is removed" << std::endl;
to_remove.push_back(nb);
}
for (std::vector<int>::const_iterator it=to_remove.begin(); it!=to_remove.end(); ++it)
else
{
int nb = *it;
point_container.erase(nb);
auto pe_uvw = list_pe_uvw.front();
std::unique_ptr<xfem::xGeomElem> geo(new xfem::xGeomElem(const_cast<AOMD::mEntity*>(pe_uvw.first)));
geo->setUVW(pe_uvw.second);
geo->SetIntegrationPointNumberForDegree(0);
geom_elem_container.insert(std::make_pair(nb, std::move(geo)));
}
}
}
for (auto id : to_remove) point_container.erase(id);
}
xExportSensors::~xExportSensors()
{
for (std::map<std::string, std::map<int, std::fstream*> >::iterator it=to_export_container.begin(); it!=to_export_container.end(); ++it)
xExportSensors::~xExportSensors()
{
for (std::map<std::string, std::map<int, std::fstream*>>::iterator it = to_export_container.begin();
it != to_export_container.end(); ++it)
{
std::map<int, std::fstream*>& filestr_map = it->second;
for (std::map<int, std::fstream*>::iterator itm = filestr_map.begin(); itm != filestr_map.end(); ++itm)
{
std::map<int, std::fstream*>& filestr_map = it->second;
for (std::map<int, std::fstream*>::iterator itm=filestr_map.begin(); itm!=filestr_map.end(); ++itm)
{
if (itm->second)
{
itm->second->close();
delete itm->second;
}
}
filestr_map.clear();
if (itm->second)
{
itm->second->close();
delete itm->second;
}
}
to_export_container.clear();
}
filestr_map.clear();
}
to_export_container.clear();
}
void xExportSensors::init(const xParseData& sensor_parser_, xfem::xMesh* mesh_, const int step_string_length_,
std::string key_label_, std::string point_label_)
{
step_string_length = step_string_length_;
std::list<std::string> string_list = sensor_parser_.getListString(key_label_);
std::list<std::string>::const_iterator it = string_list.begin();
for (; it!=string_list.end(); ++it)
{
to_export_container.insert(std::pair<std::string, std::map<int, std::fstream*> >(*it, std::map<int, std::fstream*>()));
}
std::list<xtensor::xVector<> > point_list = sensor_parser_.getListVector(point_label_);
std::list<xtensor::xVector<> >::const_iterator itp = point_list.begin();
int counter = 1;
for (; itp!=point_list.end(); ++itp, ++counter)
{
xtensor::xVector<> vec = *itp;
sensors.registerSensor(counter, Trellis_Util::mPoint(vec(0),vec(1),vec(2)), mesh_);
}
}
void xExportSensors::init(const xParseData& sensor_parser_, xfem::xMesh* mesh_, const int step_string_length_,
std::string key_label_, std::string point_label_)
{
step_string_length = step_string_length_;
std::list<std::string> string_list = sensor_parser_.getListString(key_label_);
for (auto name : string_list)
{
to_export_container.insert(std::pair<std::string, std::map<int, std::fstream*>>(name, std::map<int, std::fstream*>()));
}
int counter = 1;
for (auto pt : sensor_parser_.getListVector(point_label_)) sensors.registerSensor(counter++, pt, mesh_);
}
bool xExportSensors::toExport(const std::string& export_name_) {
std::map<std::string, std::map<int, std::fstream*> >::const_iterator it = to_export_container.find(export_name_);
if (it!=to_export_container.end()) return true;
return false;
}
bool xExportSensors::toExport(const std::string& export_name_)
{
return to_export_container.find(export_name_) != to_export_container.end();
}
} // namespace xexport
This diff is collapsed.
This diff is collapsed.
......@@ -19,7 +19,6 @@
#include "mpi.h"
// AOMD includes
#include "AOMD_SharedInfo.h"
#include "GEntity.h"
#include "mAOMD.h"
#include "mVertex.h"
......@@ -29,21 +28,21 @@
#include "xEntityFilter.h"
#include "xEntityToEntity.h"
#include "xIter.h"
#include "xOctreeGrid.h"
#include "xPartition.h"
#include "xRegularGrid.h"
#include "xUnorderedMapDataManager.h"
// xinterface aomd
#include "xAttachedDataManagerAOMD.h"
// xtensor
#include "xBoundingBox.h"
#include "xPoint.h"
struct Octree;
namespace xfem
{
class xSubMeshCreator;
class xRegularGrid;
class xRegion;
class xLevelSet;
class xSubMesh;
......@@ -59,7 +58,7 @@ typedef std::function<void(AOMD::mEntity*, xPartition&, xEntityFilter)> xGetPart
*/
AOMD::mEntity* getSource(AOMD::mEntity* e);
std::pair<xtensor::xPoint, xtensor::xPoint> compute_bounding_box(AOMD::mEntity* e);
xtensor::xBoundingBox compute_bounding_box(const AOMD::mEntity& e);
/// the xMesh adds a functionality to the mMesh class of the Aomd package.
/*!
......@@ -95,8 +94,8 @@ class xMesh //: public AOMD::mMesh
xMesh& operator=(const xMesh&) = delete;
/// return an iterator range on entity of level what.
xtool::xRange<xIter> range(int what) const { return xtool::make_range(begin(what), end(what)); }
// void modifyAllState();
// void modifyAllStateFalse();
//[[deprecated]] void modifyAllState();
//[[deprecated]] void modifyAllStateFalse();
/// remove the specific data of an xMesh relative to e. To fully Delete an entity from the xMesh, see the comment below :
void clear(AOMD::mEntity* e);
/// change the del function of the base class (call clear(e) the mMesh::del(e))
......@@ -117,6 +116,8 @@ class xMesh //: public AOMD::mMesh
void copyMesh(const xMesh& other, xinterface::aomd::xAttachedDataManagerAOMD<AOMD::mEntity*>& associated_new_entity);
/// return 2 points forming the bounding box of the mesh
xtensor::xBoundingBox compute_bounding_box() const;
//[[deprecated]]
void compute_bounding_box(xtensor::xPoint& min, xtensor::xPoint& max) const;
/// return the dimension (0, 1, 2, 3 ) of the highest level entity in the mesh (local to the proc)
......@@ -171,10 +172,20 @@ class xMesh //: public AOMD::mMesh
void cleanPartition();
void cleanClassification();
/// create an aomd octree data structure to locate elements.
/// create octree data structure to locate elements (xOctreeGrid).
void createOctree() const;
void locateElementOctree(const xtensor::xPoint& p, std::set<AOMD::mEntity*>&) const;
void locateElement(const xtensor::xPoint& p, std::set<AOMD::mEntity*>&);
/// clear the octree data structure to locate elements (xOctreeGrid).
void clearOctree() const;
/// return the octree data structure to locate elements (xOctreeGrid).
const xOctreeGrid& getOctreeGrid() const;
/// create regular grid data structure to locate elements (xRegularGrid).
void createGrid() const;
/// clear the regular grid data structure to locate elements (xRegularGrid).
void clearGrid() const;
/// return the regular grid data structure to locate elements (xRegularGrid).
const xRegularGrid& getRegularGrid() const;
std::vector<std::pair<const AOMD::mEntity*, const xtensor::xPoint>> locateElementOctree(const xtensor::xPoint& p) const;
std::vector<std::pair<const AOMD::mEntity*, const xtensor::xPoint>> locateElement(const xtensor::xPoint& p) const;
static datamanager_t<AOMD::mEntity*>& get_was_created_by();
static const datamanager_t<AOMD::mEntity*>& get_const_was_created_by();
......@@ -223,8 +234,8 @@ class xMesh //: public AOMD::mMesh
private:
partman_t part_man;
xRegularGrid* RegularGrid;
mutable std::unique_ptr<xRegularGrid> regular_grid = nullptr;
mutable std::unique_ptr<xOctreeGrid> octree_grid = nullptr;
AOMD::mMesh mesh;
// access functions to the subsets
// Subset are implemented as a map between the string that named the subset and an
......@@ -247,7 +258,7 @@ class xMesh //: public AOMD::mMesh
AOMD::mVertex* vertexInBetween(AOMD::mEntity* v0, AOMD::mEntity* v1);
void addNewTo(AOMD::mEntity* e, const std::string& side, xEntityToEntity& class_in, xEntityToEntity& class_out);
static void setDataManagers();
mutable Octree* octree;
xSubMesh* boundary;
};
////////////////////////////////// end class xMesh //////////////////////////////////////////////////
......
/*
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
and conditions.
*/
#ifndef _XOCTREEGRID_H
#define _XOCTREEGRID_H
// std
#include <forward_list>
#include <memory>
// xtensor
#include "xBoundingBox.h"
#include "xPoint.h"
#include "xVector.h"
namespace xfem
{
/// This is the template class to represent the octree search data structure.
//! The template parameters are:
//! maxElem is the maximum number of element centeroid that can be found in an octant
//! maxLevel is the maxmimum number of Level of subdivision
//! numDiv is the number of subdivision per axis. ( for a proper octree nuMDiv =2)
//! Each octant bucket of the octree represent a bounding box and entities that covers it are
//! accessible from the buckets. The octree refinment strategie is based on the number of entity
//! whose centroid are inside the bucket.
//! If this number is larger than The template parameter maxElem and the level of the bucket is
//! lower than maxLevel, the bucket is divided and the entities are moved to higher level buckets,
//! recursively, until the maxElem or maxLevel constrains are reached.
template <size_t maxElem = 5, size_t maxLevel = 20, size_t numDiv = 2>
class octantBucket
{
public:
struct elem_data
{
void *elem = nullptr; /* the pointer to a mesh Db region */
xtensor::xPoint centroid = {0., 0., 0.}; /* centroid of element bounding box inside of the octant */
xtensor::xBoundingBox bb = {{0., 0., 0.}, {0., 0., 0.}};
// \note following constructors are not need with recent compiler (gcc9) ... but it need to be spelled out at least for
// gcc4.8
elem_data(void *_elem, const xtensor::xPoint &_centroid, const xtensor::xBoundingBox &_bb)
: elem(_elem), centroid(_centroid), bb(_bb)
{
}
elem_data() = default;
};
octantBucket(const xtensor::xBoundingBox &_bb, size_t _level = 0, const octantBucket *_parent = nullptr)
: bb(_bb), level(_level), parent(_parent)
{
}
octantBucket() = default;
/// Return the level of the bucket (Number of division since root bucket)
size_t getLevel() const { return level; }
/// return the bucket Bounding box
const xtensor::xBoundingBox &getBoundingBox() const { return bb; }
/// return the list of data attached to the buckets.
const std::forward_list<elem_data> &getDatas() const { return covered_elements; }
/// call fun on all leaves of the octree, (depth first)
void visitLeaves(std::function<void(const octantBucket &buck)> fun) const;
/// return the leaf bucket containing p if exist it exist, nullptr otherwise
const octantBucket<maxElem, maxLevel, numDiv> *findBucket(const xtensor::xPoint &p) const;
/// Add element to leaves octant bucket's list if edata.bb cover the leaf boundingbox.
//! if a leaf has more than maxElem element centeroid, the leaf is subdivided and all data are moved to these new leaves
//! recursivelly
void addElement(const elem_data &edata);
private:
xtensor::xBoundingBox bb; // bounding box.
size_t level = 0; // level in the octree of the bucket
const octantBucket *parent = nullptr; // link to the parent bucket */
// number of element centeroid contains in the bucket.
size_t contained_centroids_size = 0;
// list of element whose BB intersect octant.
std::forward_list<elem_data> covered_elements;
static constexpr size_t numChilds = numDiv * numDiv * numDiv;
// Pointer to all the childs (numBucks ... ) if any.
std::unique_ptr<std::array<octantBucket, numChilds>> childs;
/// try to refine
//! prerequist : bucket has no child
//! Returns true for success, false for failure (no memory left).
bool subdivideOctantBucket();
};
using xOctreeGrid = octantBucket<5, 20, 2>;
template <size_t maxElem, size_t maxLevel, size_t numDiv>
void octantBucket<maxElem, maxLevel, numDiv>::visitLeaves(
std::function<void(const octantBucket<maxElem, maxLevel, numDiv> &buck)> fun) const
{
if (childs)
{
for (const auto &child : *childs) child.visitLeaves(fun);
return;
}
fun(*this);
}
template <size_t maxElem, size_t maxLevel, size_t numDiv>
const octantBucket<maxElem, maxLevel, numDiv> *octantBucket<maxElem, maxLevel, numDiv>::findBucket(const xtensor::xPoint &p) const
{
if (!bb.contains(p)) return nullptr;
if (!childs) return this;
for (const octantBucket &child : *childs)
{
const octantBucket *found = child.findBucket(p);
if (found) return found;
}
return nullptr;
}
template <size_t maxElem, size_t maxLevel, size_t numDiv>
void octantBucket<maxElem, maxLevel, numDiv>::addElement(const elem_data &edata)
{
if (!edata.bb.cover(bb)) return;
if (childs)
{
for (octantBucket &bucket : *childs) bucket.addElement(edata);
return;
}
// At this point, I'am on a leaf that cover the bounding box of elem.
// checking if edata already in the list. Do nothing if its the case.
for (const auto &data : covered_elements)
{
if (data.elem == edata.elem)
{
if (!(data.bb == edata.bb && data.centroid == edata.centroid))
std::cout << "Warning !!! added same element with diff centeroid and bb" << std::endl;
return;
}
}
covered_elements.push_front(edata);
if (bb.contains(edata.centroid)) ++contained_centroids_size;
if (contained_centroids_size <= maxElem) return;
if (level < maxLevel)
subdivideOctantBucket();
else
std::cout << "Warning : maxLevel Reach in xOctreeGrid, more than " << maxElem << " In currrent octant " << __FILE__
<< " line " << __LINE__ << std::endl;
return;
}
template <size_t maxElem, size_t maxLevel, size_t numDiv>
bool octantBucket<maxElem, maxLevel, numDiv>::subdivideOctantBucket()
{
if (childs) throw;
// std::make_unique is c++14
// childs = std::make_unique<std::array<octantBucket, numChilds>>();
childs = std::unique_ptr<std::array<octantBucket, numChilds>>(new std::array<octantBucket, numChilds>);
if (!childs)
{
std::cerr << "Error, subdivideOctantBucket could not allocate enough space" << std::endl;
return false;
}
const xtensor::xPoint tmp = (bb.max - bb.min) / numDiv;
for (size_t k = 0; k < numDiv; ++k)
{
size_t shift = k * numDiv * numDiv;
for (int j = 0; j < numDiv; ++j)
{
for (int i = 0; i < numDiv; ++i)
{
octantBucket &child = (*childs)[shift + i];
child.parent = this;
child.level = level + 1;
child.bb = {bb.min + xtensor::xPoint(tmp[0] * i, tmp[1] * j, tmp[2] * k),
bb.min + xtensor::xPoint(tmp[0] * (i + 1), tmp[1] * (j + 1), tmp[2] * (k + 1))};
}
shift += numDiv;
}
}
contained_centroids_size = 0;
for (auto &child : *childs)
for (const elem_data &data : covered_elements) child.addElement(data);
covered_elements.clear();
return true;
}
} // namespace xfem
#endif
......@@ -15,177 +15,118 @@ namespace xfem
{
using AOMD::mEntity;
using AOMD::mVertex;
using xtensor::xBoundingBox;
using xtensor::xPoint;
xRegularGrid::xRegularGrid(xMesh *m, int nx, int ny, int nz)
xRegularGrid::xRegularGrid(const xtensor::xBoundingBox &bb, int nx, int ny, int nz) : N{nx, ny, nz}, numbox{nx * ny * nz}
{
N[0] = nx;
N[1] = ny;
N[2] = nz;
numbox = nx * ny * nz;
xtensor::xPoint mini, maxi;
m->compute_bounding_box(mini, maxi);
for (int i = 0; i < 3; i++)
{
Pmin(i) = mini(i) - TOL;
Pmax(i) = maxi(i) + TOL;
}
for (int i = 0; i < numbox; i++)
{
xBrick b(i);
Table.push_back(b);
}
initTable();
const xtensor::xPoint PTOL{TOL, TOL, TOL};
Pmin = bb.min - PTOL;
Pmax = bb.max + PTOL;
}
for (AOMD::mEntity *e : m->range(m->dim()))
xRegularGrid::xRegularGrid(const xMesh &m, int nx, int ny, int nz) : N{nx, ny, nz}, numbox{nx * ny * nz}
{
initTable();
const xBoundingBox mbb = m.compute_bounding_box();
const xtensor::xPoint PTOL{TOL, TOL, TOL};
Pmin = mbb.min - PTOL;
Pmax = mbb.max + PTOL;
const auto &mappingbuilder = xMappingBuilderHolderSingleton::instance();
for (const AOMD::mEntity *e : m.range(m.dim()))
{
xtensor::xPoint min, max;
BoundingBoxEntity(e, min, max);
AddObject(e, min, max);
std::unique_ptr<xmapping::xMapping> pmapping(mappingbuilder.buildMapping(*e));
addObject(*e, pmapping->boundingBox());
}
}
void xRegularGrid ::Print()
void xRegularGrid::print(std::ostream &out) const
{
printf("Here are the infos in the regular grid\n");
printf("#####################################\n");
int index, nb, i, j, k, l;
xBrick *pBrick;
mEntity *e;
for (i = 0; i < N[0]; i++)
{
for (j = 0; j < N[1]; j++)
{
for (k = 0; k < N[2]; k++)
out << "Here are the infos in the regular grid\n";
out << "#####################################\n";
for (int i = 0; i < N[0]; ++i)
for (int j = 0; j < N[1]; ++j)
for (int k = 0; k < N[2]; ++k)
{
index = i + j * N[0] + k * N[0] * N[1];
pBrick = getBrick(index);
nb = pBrick->Objects.size();
printf("The brick at location %d, %d, %d with index %d contains %d elts\n", i, j, k, index, nb);
int index = getIndex(i, j, k);
const xBrick &brick = getBrick(index);
int nb = brick.Objects.size();
out << "The brick at location " << i << " " << j << " " << k << " with index " << index << " contains " << nb
<< "elts\n";
if (nb)
{
printf("The list is ");
for (l = 0; l < nb; l++)
{
e = pBrick->Objects[l];
printf("%d ", e->getId());
}
printf("\n");
out << "The list is ";
for (const AOMD::mEntity *e : brick.Objects) out << e->getId() << " ";
out << "\n";
}
}
}
}
printf("End of printing the regular grid infos\n");
printf("#####################################\n");
out << "End of printing the regular grid infos\n";
out << "#####################################\n";
return;
}