* @param numPts number of vertices
*
*/
- BoundingBox::BoundingBox(const double** pts, const int numPts);
+ BoundingBox::BoundingBox(const double** pts, const int numPts)
{
using namespace std;
assert(numPts > 1);
* Constructor creating box from union of two boxes
*
*/
- BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box1)
+ BoundingBox::BoundingBox(const BoundingBox& box1, const BoundingBox& box2)
{
+ using namespace std;
for(BoxCoord c = XMIN ; c <= ZMIN ; c = BoxCoord(c + 1))
{
- _coords[c] = min(box1[c], box2[c]);
- _coords[c + 3] = max(box1[c + 3], box2[c + 3]);
+ _coords[c] = min(box1._coords[c], box2._coords[c]);
+ _coords[c + 3] = max(box1._coords[c + 3], box2._coords[c + 3]);
}
}
* Constructor creating box from union of two boxes
*
*/
- BoundingBox(const BoundingBox& box1, const BoundingBox& box1) {}
+ BoundingBox(const BoundingBox& box1, const BoundingBox& box2);
/**
* Destructor
void updateWithPoint(const double* pt);
private:
- double[6] _coords;
+ double _coords[6];
};
}
else // recursion
{
- typedef MeshElement::ZoneBBoxComparison Cmp;
+
// std::cout << " - Recursion" << std::endl;
RegionNode* leftNode = new RegionNode();
// split current source region
//} decide on axis
- static Cmp::BoxAxis axis = Cmp::X;
+ static BoundingBox::BoxCoord axis = BoundingBox::XMAX;
currNode->getSrcRegion().split(leftNode->getSrcRegion(), rightNode->getSrcRegion(), axis);
// ugly hack to avoid problem with enum which does not start at 0
// I guess I ought to implement ++ for it instead ...
// Anyway, it basically chooses the next axis, circually
- axis = (axis != Cmp::Z) ? static_cast<Cmp::BoxAxis>(axis + 1) : Cmp::X;
+ axis = (axis != BoundingBox::ZMAX) ? static_cast<BoundingBox::BoxCoord>(axis + 1) : BoundingBox::XMAX;
// add target elements of curr node that overlap the two new nodes
// std::cout << " -- Adding target elements" << std::endl;
#include "MeshElement.hxx"
+#include "TetraAffineTransform.hxx"
+#include "TransformedTriangle.hxx"
+
namespace INTERP_UTILS
{
* @param type geometric type of the element
* @param index global number of element in the mesh
*/
- MeshElement::MeshElement(const int index, const MEDEN::medGeometryElement type, const MEDMEM::Mesh& mesh)
+ MeshElement::MeshElement(const int index, const MED_EN::medGeometryElement type, const MEDMEM::MESH& mesh)
: _index(index - 1), _mesh(&mesh), _type(type), _box(0)
{
// get coordinates of vertices
const double* MeshElement::getCoordsOfNode(int node) const
{
const int nodeOffset = node - 1;
- const int elemIdx = mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index];
+ const int elemIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_CELL)[_index];
return &(_mesh->getCoordinates(MED_FULL_INTERLACE)[elemIdx + nodeOffset]);
}
*/
void MeshElement::triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const
{
- // not implemented
+ switch(_type)
+ {
+ case MED_TETRA4 :
+ // triangles == faces
+ const int beginFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index];
+ const int endFaceIdx = _mesh->getConnectivityIndex(MED_DESCENDING, MED_CELL)[_index + 1];
+
+ for(int i = beginFaceIdx ; i < endFaceIdx ; ++i) // loop over faces of element
+ {
+ const int faceIdx = _mesh->getConnectivity(MED_FULL_INTERLACE, MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS)[i - 1];
+ const int beginNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx - 1];
+ const int endNodeIdx = _mesh->getConnectivityIndex(MED_NODAL, MED_FACE)[faceIdx];
+
+ double transformedPts[9];
+
+ assert(endNodeIdx - beginNodeIdx == 3);
+
+ for(int j = 0 ; j < 3 ; ++j) // loop over nodes of face
+ {
+ // { optimise here using getCoordinatesForNode ?
+ // could maybe use the connNodeIdx directly and only transform each node once
+ // instead of once per face
+ const int connNodeIdx =
+ _mesh->getConnectivity(MED_FULL_INTERLACE, MED_NODAL, MED_FACE, MED_ALL_ELEMENTS)[beginNodeIdx + j - 1];
+ const double* pt = &(_mesh->getCoordinates(MED_FULL_INTERLACE)[connNodeIdx - 1]);
+
+ // transform
+ T.apply(&transformedPts[3*j], pt);
+ }
+
+ triangles.push_back(TransformedTriangle(&transformedPts[0], &transformedPts[3], &transformedPts[6]));
+ }
+
+ break;
+ default:
+ break;
+ }
+ }
+
+ int MeshElement::getIndex() const
+ {
+ return _index + 1;
+ }
+
+
+ /// ElementBBoxOrder
+ bool ElementBBoxOrder::operator()( MeshElement* elem1, MeshElement* elem2)
+ {
+ assert(elem1 != 0);
+ assert(elem2 != 0);
+ assert(elem1->_box != 0);
+ assert(elem2->_box != 0);
+
+ const double coord1 = elem1->_box->getCoordinate(_coord);
+ const double coord2 = elem2->_box->getCoordinate(_coord);
+
+ return coord1 < coord2;
}
+
};
#ifndef __MESH_ELEMENT_HXX__
#define __MESH_ELEMENT_HXX__
+#include "MEDMEM_define.hxx"
+#include "MEDMEM_Mesh.hxx"
+
+#include "BoundingBox.hxx"
+
+using namespace MEDMEM;
+using namespace MED_EN;
+
+
+
+
namespace INTERP_UTILS
{
+ class TransformedTriangle;
+ class TetraAffineTransform;
+
/**
* Class representing a single element of a mesh together with its bounding box.
* It permits access to the nodes of the element, to test for disunion and inclusion with other elements
*/
class MeshElement
{
- public:
- /**
- * Constructor
- *
- * @param mesh mesh that the element belongs to
- * @param type geometric type of the element
- * @param index connectivity index of element in the mesh
- */
- MeshElement(const int index, const MEDEN::medGeometryElement type, const MEDMEM::Mesh& mesh);
-
- /**
- * Destructor
- *
- */
- ~MeshElement();
-
- /**
- * Determines if this element is in the interior of another element
- * by calculating the triple products for each point of this element with respect
- * to all the faces of the other object (faces must be triangulated ... ) If all triple
- * products have the same sign, then the element is in the interior of the other element
- *
- * @param otherElement the supposedly enclosing element
- * @returns true if this element is enclosed in the other element, false if not
- */
- bool isElementIncludedIn(const MeshElement& otherElement) const;
-
- /**
- * Determines whether the intersection of this element is trivially empty. This is done by checking for each
- * face of one element if it is such that all the vertices of the other element is on the same side of this face.
- * If there is such a face, then the intersection is trivially empty. If there is no such face, we do not know if
- * the intersection is empty.
- *
- * @pre The elements are convex. If this is no true, we return false.
- * @param otherElement the element believed to be disjoint with this one
- * @returns true if the two elements are convex and there exists a face of this element such as described
- * above, false otherwise
- */
- bool isElementTriviallyDisjointWith(const MeshElement& otherElement) const;
- /**
- * Returns the number of nodes of this element
- *
- * @returns the number of nodes of this element
- */
- int getNumberNodes() const;
-
- /**
- * Returns the coordinates of a node of this element
- * (1 <= node <= #nodes)
- *
- * @param node the node for which the coordinates are sought
- * @returns pointer to an array of 3 doubles containing the coordinates
- */
- const double* getCoordsOfNode(int node);
+ friend class ElementBBoxOrder;
+ friend class MeshRegion;
+ public:
+
/**
* Constructor
*
* @param type geometric type of the element
* @param index connectivity index of element in the mesh
*/
- MeshElement(const int index, const MEDEN::medGeometryElement type, const MEDMEM::Mesh& mesh);
+ MeshElement(const int index, const MED_EN::medGeometryElement type, const MEDMEM::MESH& mesh);
/**
* Destructor
* @param node the node for which the coordinates are sought
* @returns pointer to an array of 3 doubles containing the coordinates
*/
- const double* getCoordsOfNode(int node);
+ const double* getCoordsOfNode(int node) const;
/**
* Triangulate the faces of this element and apply an affine Transform to the triangles
* @param triangles vector in which triangles are stored
* @param T affine transform that is applied to the nodes of the triangles
*/
- void triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T);
+ void triangulate(std::vector<TransformedTriangle>& triangles, const TetraAffineTransform& T) const;
+
+ int getIndex() const;
private:
const int _index;
- const MEDMEM::Mesh* _mesh;
- const MEDEN::medGeometryElement _type;
- BoundingBox _box; // should not change after initialisation
+ const MEDMEM::MESH* _mesh;
+ const MED_EN::medGeometryElement _type;
+ BoundingBox* _box; // should not change after initialisation
};
- // to be added : triangulation methods
-};
+ class ElementBBoxOrder
+ {
+ public :
+
+ ElementBBoxOrder(BoundingBox::BoxCoord coord)
+ : _coord(coord)
+ {
+ }
+
+ bool operator()(MeshElement* elem1, MeshElement* elem2);
+
+ private :
+ BoundingBox::BoxCoord _coord;
+ };
+};
#endif
#include "MeshRegion.hxx"
+#include "MeshElement.hxx"
+
+//#include <math.h>
+
namespace INTERP_UTILS
{
+
/**
- * Class representing a set of elements in a mesh together with their bounding box.
- * It permits to split itself in two, which is used in the depth-first search filtering process.
- *
+ * Default constructor
+ *
*/
- class MeshRegion
+ MeshRegion::MeshRegion()
+ : _box(0)
{
- public:
-
- /**
- * Default constructor
- *
- */
- MeshRegion::MeshRegion()
- : _box(0)
- {
- }
+ }
- /**
- * Destructor
- *
- */
- MeshRegion::~MeshRegion()
- {
- if(_box != 0)
- {
- delete _box;
- }
-
- /**
- * Adds an element to the region, updating the bounding box.
- *
- * @param element pointer to element to add to region
- *
- */
- void MeshRegion::addElement(MeshElement* const element)
+ /**
+ * Destructor
+ *
+ */
+ MeshRegion::~MeshRegion()
+ {
+ if(_box != 0)
{
- _elements.push_back(element);
+ delete _box;
+ }
+ }
+
+ /**
+ * Adds an element to the region, updating the bounding box.
+ *
+ * @param element pointer to element to add to region
+ *
+ */
+ void MeshRegion::addElement(MeshElement* const element)
+ {
+ _elements.push_back(element);
- if(_box == 0)
+ if(_box == 0)
+ {
+ const int numNodes = element->getNumberNodes();
+ const double* pts[numNodes];
+
+ // get coordinates of elements
+ for(int i = 0 ; i < numNodes ; ++i)
{
- // get coordinates of elements
- _box = new BoundingBox();
+ pts[i] = element->getCoordsOfNode(i + 1);
}
+
+ _box = new BoundingBox(pts, numNodes);
+
+ }
+ }
+
+ /**
+ * Splits the region in two along the given axis, copying the elements with bounding boxes whose maximum
+ * coordinate along the axis are smaller than the middle of the bounding box of this region in region1. The
+ * rest of the elements are copied to region2.
+ *
+ * @param region1 region in which to store one half of this region
+ * @param region2 region in which to store the other of this region
+ * @param coord coordinate of BoundingBox to use when splitting the region
+ *
+ */
+ void MeshRegion::split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord)
+ {
+ ElementBBoxOrder cmp(coord);
+
+ // sort elements by their bounding boxes
+ std::sort(_elements.begin(), _elements.end(), cmp);
+
+ // put the first half of the elements in region1 and the
+ // rest in region2
+ std::vector<MeshElement*>::const_iterator iter = _elements.begin();
+ int elemCount = 0;
+
+ while(elemCount < static_cast<int>(_elements.size() / 2))
+ {
+ region1.addElement(*iter);
+ ++iter;
+ ++elemCount;
+ }
+
+ while(iter != _elements.end())
+ {
+ region2.addElement(*iter);
+ ++iter;
}
- /**
- * Splits the region in two along the given axis, copying the elements with bounding boxes whose maximum
- * coordinate along the axis are smaller than the middle of the bounding box of this region in region1. The
- * rest of the elements are copied to region2.
- *
- * @param region1 region in which to store one half of this region
- * @param region2 region in which to store the other of this region
- * @param axis axis along which to split the region
- *
- */
- void MeshRegion::split(Region& region1, Region& region2, int axis) const;
-
- private:
- // Vector of pointers to elements. NB : these pointers are not owned by the region object, and are thus
- // neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class
- std::vector<MeshElement*> _elements;
- BoundingBox* _box;
-
- };
+ // assert(std::abs(region1._elements.size() - region2._elements.size()) < 2);
+
+ }
+
+ bool MeshRegion::isDisjointWithElementBoundingBox(const MeshElement& elem) const
+ {
+ assert(_box != 0);
+ assert(elem._box != 0);
+
+ return _box->isDisjointWith(*(elem._box));
+ }
+ std::vector<MeshElement*>::const_iterator MeshRegion::getBeginElements() const
+ {
+ return _elements.begin();
+ }
+
+ std::vector<MeshElement*>::const_iterator MeshRegion::getEndElements() const
+ {
+ return _elements.end();
+ }
+
+ int MeshRegion::getNumberOfElements() const
+ {
+ return _elements.size();
+ }
};
-#endif
+
#include <vector>
+#include "BoundingBox.hxx"
+
namespace INTERP_UTILS
{
+ class MeshElement;
+
/**
* Class representing a set of elements in a mesh together with their bounding box.
* It permits to split itself in two, which is used in the depth-first search filtering process.
* @param axis axis along which to split the region
*
*/
- void split(Region& region1, Region& region2, int axis) const;
+ void split(MeshRegion& region1, MeshRegion& region2, BoundingBox::BoxCoord coord);
+
+ bool isDisjointWithElementBoundingBox(const MeshElement& elem) const;
+
+ std::vector<MeshElement*>::const_iterator getBeginElements() const;
+
+ std::vector<MeshElement*>::const_iterator getEndElements() const;
+
+ int getNumberOfElements() const;
private:
// Vector of pointers to elements. NB : these pointers are not owned by the region object, and are thus
// neither allocated or liberated in this class. The elements must therefore be allocated and liberated outside this class
std::vector<MeshElement*> _elements;
- BoundingBox _box;
+ BoundingBox* _box;
};
--- /dev/null
+#include "RegionNode.hxx"
+
+namespace INTERP_UTILS
+{
+ /**
+ * Default constructor
+ *
+ */
+ RegionNode::RegionNode()
+ {
+ }
+
+ /**
+ * Destructor
+ *
+ */
+ RegionNode::~RegionNode()
+ {
+ }
+
+
+ /**
+ * Accessor to source region
+ *
+ * @return reference to source region
+ */
+ MeshRegion& RegionNode::getSrcRegion()
+ {
+ return _srcRegion;
+ }
+ /**
+ * Accessor to target region
+ *
+ * @return reference to target region
+ */
+ MeshRegion& RegionNode::getTargetRegion()
+ {
+ return _targetRegion;
+ }
+
+};
#ifndef __REGION_NODE_HXX__
#define __REGION_NODE_HXX__
+#include "MeshRegion.hxx"
+
namespace INTERP_UTILS
{
+
/**
* Class containing a tuplet of a source region and a target region.
* This is used as the object to put on the stack in the depth-first search
~RegionNode();
/**
+ * Accessor to source region
*
+ * @return reference to source region
*/
- MeshRegion& getSrcRegion() const;
+ MeshRegion& getSrcRegion();
/**
+ * Accessor to target region
*
+ * @return reference to target region
*/
- MeshRegion& getTargetRegion() const;
+ MeshRegion& getTargetRegion();
private:
