Geometry2D Class Reference

#include <geometry2d.h>

Inheritance diagram for Geometry2D:

Collaboration diagram for Geometry2D:

Public Member Functions
	Geometry2D (const MInt solverId_, const MPI_Comm comm)
	Geometry2D (const MInt solverId_, const MString filename, const MPI_Comm comm)
	~Geometry2D ()
virtual MInt	getLineIntersectionElements (MFloat *targetRegion, std::vector< MInt > &nodeList)
	Returns the ids of all elements, cut by a orthogonal line (or rectangular region) More...
virtual MInt	getLineIntersectionElementsOld1 (MFloat *targetRegion, std::vector< MInt > &nodeList)
virtual MInt	getLineIntersectionElementsOld2 (MFloat *targetRegion, MInt *spaceDirection, std::vector< MInt > &nodeList)
virtual MInt	getIntersectionElements (MFloat *targetRegion, std::vector< MInt > &nodeList)
	Determines all elements that are inside or intersect the target region. More...
virtual MInt	getIntersectionElements (MFloat *targetRegion, std::vector< MInt > &nodeList, MFloat cellHalfLength, const MFloat *const cell_coords)
	Detemines whether a triangle and a gridcell intersect with the separating axis theorem (SAT) More...
virtual MBool	edgeTriangleIntersection (MFloat *trianglePoint1, MFloat *trianglePoint2, MFloat *trianglePoint3, MFloat *edgePoint1, MFloat *edgePoint2)
	Determine intersection between an edge and a triangle (in 2D between 2 edges) More...
virtual MInt	getIntersectionMBElements (MFloat *targetRegion, std::vector< MInt > &nodeList)
virtual MInt	getLineIntersectionMBElements (MFloat *targetRegion, std::vector< MInt > &nodeList)
virtual MInt	getSphereIntersectionMBElements (MFloat *P, MFloat radius, std::vector< MInt > &nodeList)
virtual void	MoveAllMBElementVertex (MFloat *dx)
virtual void	MoveMBElementVertex (MInt e, MInt v, MFloat *dx)
virtual void	ReplaceMBElementVertex (MInt e, MInt v, MFloat *np)
virtual void	UpdateMBNormalVector (MInt)
virtual void	UpdateMBBoundingBox ()
virtual void	UpdateADT ()
Public Member Functions inherited from Geometry< 2 >
	Geometry (const MInt solverId_, const MPI_Comm comm)
virtual	~Geometry ()=default
MInt	solverId () const
MPI_Comm	mpiComm () const
MInt	domainId () const
MInt	noDomains () const
GeometryContext &	geometryContext ()
void	getBoundingBox (MFloat *const bBox) const
	Returns the bounding box for the geometry. More...
const MFloat *	boundingBox () const
void	setHaloElementOffset (MInt off)
MInt	getHaloElementOffset () const
MBool	isOnGeometry (const MFloat, const MFloat *, MString)
MBool	vectorsEqual (MFloat *a, MFloat *b)
	Compares two vectors entry by entry. More...
MFloat	calcCircumference (MFloat **bndVs, MInt num)
	Returns the circumference of a segment. More...
virtual MInt	boundaryCheck (MFloat *, MFloat, MFloat *, MInt *)
virtual MInt	getIntersectionElements (MFloat *, std::vector< MInt > &)
virtual MInt	getIntersectionElements (MFloat *, std::vector< MInt > &, MFloat, const MFloat *const)
virtual MInt	getLineIntersectionElementsOld1 (MFloat *, std::vector< MInt > &)
virtual MInt	getLineIntersectionElementsOld2 (MFloat *, MInt *, std::vector< MInt > &)
virtual MInt	getLineIntersectionElements (MFloat *, std::vector< MInt > &)
virtual MInt	getLineIntersectionElements (MFloat *)
virtual MBool	getClosestLineIntersectionLength (MInt, const std::vector< MInt > &, MFloat *, MFloat *)
void	getLineIntersectingElementsBcIds (const MFloat *const line, std::set< MInt > &bcIds)
	Return the set of boundary condition ids of the elements cut by the given line. More...
virtual MBool	edgeTriangleIntersection (MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
virtual MBool	edgeTriangleIntersectionLB (MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
virtual MBool	getLineTriangleIntersectionSimple (MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
virtual MBool	getLineTriangleIntersectionSimpleDistance (MFloat *, MFloat *, MFloat *, MFloat *, MFloat *, MFloat *)
virtual MBool	getLineTriangleIntersection (const MFloat *const, const MFloat *const, const MFloat, const MFloat *const, const MFloat *const, const MFloat *const, MFloat *, MFloat *, MFloat *, MFloat *)
void	getBoundingBoxMB (MFloat *const bBox) const
virtual MInt	getIntersectionMBElements (MFloat *, std::vector< MInt > &)
virtual MInt	getLineIntersectionMBElements (MFloat *, std::vector< MInt > &)
virtual MInt	getLineIntersectionMBElements2 (MFloat *, MInt *, std::vector< MInt > &, MInt)
virtual MInt	getSphereIntersectionMBElements (MFloat *, MFloat, std::vector< MInt > &)
virtual void	MoveAllMBElementVertex (MFloat *)
virtual void	MoveMBElementVertex (MInt, MInt, MFloat *)
virtual void	ReplaceMBElementVertex (MInt, MInt, MFloat *)
virtual void	UpdateMBNormalVector (MInt)
virtual void	UpdateMBBoundingBox ()
virtual void	UpdateADT ()
virtual void	collectGlobalMemoryUsage ()
virtual void	writeSTL (const MChar *)
virtual void	writeADTAndSTLToNetCDF (const MChar *)
virtual void	writeSTLMB (const MChar *, MInt &, MInt *&)
virtual void	readSTLNetCDF (const MChar *)
virtual void	logStatistics ()
virtual MInt	GetNoElements ()
virtual MInt	GetNoSegments ()
virtual MInt *	GetBoundaryIds (MInt *noAllBcs)
virtual MFloat **	GetBoundaryVertices (MInt, MFloat *, MInt *, MInt, MInt *)
virtual MBool	isEdgeAlreadyInCollection (std::vector< std::pair< MFloat *, MFloat * > >, MFloat *, MFloat *, MInt *)
virtual MFloat	GetBoundarySize (MInt)
virtual MFloat	GetBoundarySize (MFloat *, MInt *, MInt)
virtual void	determineSegmentOwnership (MInt, MInt *, MInt *, MInt *, MInt *)
virtual MFloat	getBndMaxRadius (MFloat **, MInt)
virtual void	rebuildAdtTree ()
virtual void	calculateBoundingBox ()
virtual void	writeParallelGeometryVTK (MString)
virtual void	addElement (MFloat *)
virtual void	copyElement (MInt, MInt)
virtual void	resizeCollector (MInt)
MBool	pointIsInside (const MFloat *const coordinates)
	Determines if a point is inside of the geometry. More...
MBool	pointIsInside (const MFloat *const coordinates, MInt *numcutsperdir)
	Determines if a point is inside of the geometry. More...
MBool	pointIsInside2 (const MFloat *const coordinates, MInt *numcutsperdir=nullptr)
	Determines if a point is in or outside the geometry. More...
MBool	pointIsInsideMBElements (const MFloat *const coordinates, MInt *, MInt *, MInt)
MBool	pointIsInsideMBElements2 (const MFloat *const coordinates, MInt *, MInt *, MInt)
void	determineRayIntersectedElements (const MFloat *const coordinates, std::vector< std::vector< MInt > > *resultnodes)
	returns the geometry elements that have a cut with rays originating in the provided coordinates More...
MInt	noBoundaryIds ()

Protected Member Functions
virtual void	readSegments ()
void	countSegmentLinesASCII (const MString &fileName, MInt *noElements)
	counts the number of lines in an ASCII file, should be the in the first line More...
void	readSegmentLinesASCII (MString fileName, Collector< element< 2 > > *elemCollector, MInt bndCndId, MInt segmentId, MInt *offset)
	reads the lines in an ASCII file More...
void	calculateBoundingBox ()
	Calculates the global bounding box of the geometry. More...
Protected Member Functions inherited from Geometry< 2 >
virtual void	readSegments ()
virtual void	writeSegmentsToDX ()

Protected Attributes
MBool	m_GFieldInitFromSTL
MInt	m_levelSetIntfBndId {}
MInt *	m_levelSetIntfBndIds {}
MInt	m_noLevelSetIntfBndIds {}
Protected Attributes inherited from Geometry< 2 >
MInt	m_noSegments
std::array< MFloat, 2 *nDim >	m_minMax
MInt	m_noElements
MString	m_segmentBaseName
bodyMap	m_bodyMap
bodyIterator	m_bodyIt
MInt	m_noMBElements
MInt	m_noBoundaryIds
MInt *	m_boundaryIds
MInt *	m_allBCs
MInt	m_noAllBCs
MBool	m_flowSolver

Static Protected Attributes
static constexpr const MInt	nDim = 2

Additional Inherited Members
Public Attributes inherited from Geometry< 2 >
std::vector< MInt >	m_segmentOffsets
std::vector< MInt >	m_segmentOffsetsWithoutMB
Collector< element< nDim > > *	m_elements
element< nDim > *	elements
GeometryAdt< nDim > *	m_adt
Collector< element< nDim > > *	m_mbelements
element< nDim > *	mbelements
std::array< MFloat, 2 *nDim >	m_mbminMax
MFloat	m_mbMidPnt [3]
MBool *	m_ownSegmentId
std::set< MInt >	m_uniqueOriginalTriId
MFloat	m_parGeomMemFactor
MString	m_inOutTest
MBool	m_parallelGeometry
MBool	m_debugParGeom
MString	m_parallelGeomFileName

Detailed Description

Specialized 2D implementation

Definition at line 15 of file geometry2d.h.

Constructor & Destructor Documentation

Geometry2D() [1/2]

Geometry2D::Geometry2D	(	const MInt	solverId_,
		const MPI_Comm	comm
	)

Definition at line 17 of file geometry2d.cpp.

                                                                : Geometry<2>(solverId_, comm) {
  TRACE();
 
  m_adt = 0;
  m_elements = 0;
  m_noElements = 0;
 
  // moving boundary
  m_mbelements = 0;
  m_noMBElements = 0;
  m_GFieldInitFromSTL = 0;
  m_levelSetIntfBndId = 0;
  m_noLevelSetIntfBndIds = 0;
  m_GFieldInitFromSTL = Context::getSolverProperty<MBool>("GFieldInitFromSTL", solverId(), AT_, &m_GFieldInitFromSTL);
 
  if(m_GFieldInitFromSTL) {
    m_levelSetIntfBndId = Context::getSolverProperty<MInt>("levelSetIntfBndId", solverId(), AT_, &m_levelSetIntfBndId);
    if(Context::propertyExists("bodyBndryCndIds", solverId())
       || Context::propertyExists("GFieldInitFromSTLBndCndIds", solverId())) {
      const MInt noBodyBndIds = Context::propertyLength("bodyBndryCndIds", solverId());
      const MInt noBndIds = Context::propertyLength("GFieldInitFromSTLBndCndIds", solverId());
 
      m_noLevelSetIntfBndIds = noBodyBndIds + noBndIds;
      m_levelSetIntfBndIds = new MInt[m_noLevelSetIntfBndIds];
 
      for(MInt i = 0; i < noBodyBndIds; i++) {
        m_levelSetIntfBndIds[i] = Context::getSolverProperty<MInt>("bodyBndryCndIds", solverId(), AT_, i);
        if(domainId() == 0) {
          cerr << " levelSetIntfBndIds: " << m_levelSetIntfBndIds[i] << endl;
        }
      }
 
      for(MInt i = noBodyBndIds; i < m_noLevelSetIntfBndIds; i++) {
        m_levelSetIntfBndIds[i] =
            Context::getSolverProperty<MInt>("GFieldInitFromSTLBndCndIds", solverId(), AT_, i - noBodyBndIds);
        if(domainId() == 0) {
          cerr << " levelSetIntfBndIds: " << m_levelSetIntfBndIds[i] << endl;
        }
      }
 
      if(domainId() == 0) {
        cerr << " noLevelSetIntfBndIds: " << m_noLevelSetIntfBndIds << endl;
      }
 
    } else {
      m_levelSetIntfBndIds = nullptr;
      //       for(MInt i=0; i < m_noLevelSetIntfBndIds; i++){
      //         m_levelSetIntfBndIds[i] = m_levelSetIntfBndId;
      //       }
    }
  }
 
  MString testcaseDir = "./";
  testcaseDir = Context::getSolverProperty<MString>("testcaseDir", solverId(), AT_, &testcaseDir);
  MString inputDir = "./";
  inputDir = Context::getSolverProperty<MString>("inputDir", solverId(), AT_, &inputDir);
 
  MString tmpFileName =
      testcaseDir + inputDir + Context::getSolverProperty<MString>("geometryInputFileName", solverId(), AT_);
 
  if(tmpFileName.find(".toml") == MString::npos) {
    geometryContext().readPropertyFile(NETCDF, tmpFileName.c_str());
  } else {
    geometryContext().readPropertyFile(TOML, tmpFileName.c_str());
  }
  m_bodyMap = geometryContext().getBodies();
  readSegments();
  m_adt = new GeometryAdt<2>(this);
  m_adt->buildTree();
 
  // moving boundary
  if(m_GFieldInitFromSTL) {
    m_adt->buildTreeMB();
  }
}

Geometry2D() [2/2]

Geometry2D::Geometry2D	(	const MInt	solverId_,
		const MString	filename,
		const MPI_Comm	comm
	)

Definition at line 98 of file geometry2d.cpp.

                                                                                        : Geometry(solverId_, comm) {
  TRACE();
 
  m_log << "reading the auxiliary stl " << filename << endl;
 
  m_GFieldInitFromSTL = 0;
 
  m_noElements = 0;
  m_noMBElements = 0;
 
  countSegmentLinesASCII(filename, &m_noElements);
 
  m_log << "  number of Elements: " << m_noElements << endl;
 
  m_elements = new Collector<element<nDim>>(m_noElements, nDim, 0);
 
  MInt bla = 0;
  readSegmentLinesASCII(filename, m_elements, 0, 0, &bla);
 
  elements = &(m_elements->a[0]);
 
  Geometry2D::calculateBoundingBox();
 
  m_adt = new GeometryAdt<2>(this);
  m_adt->buildTree();
}

~Geometry2D()

Geometry2D::~Geometry2D

(

)

Definition at line 129 of file geometry2d.cpp.

                        {
  TRACE();
 
  delete m_adt;
  delete m_elements;
 
  if(m_mbelements) delete m_mbelements;
}

Member Function Documentation

calculateBoundingBox()

void Geometry2D::calculateBoundingBox ( )

protectedvirtual

Author

Thomas Schilden

Date

27.12.2016

Reimplemented from Geometry< 2 >.

Definition at line 1108 of file geometry2d.cpp.

                                      {
  TRACE();
 
  // Calculate bounding box of segment
  for(MInt j = 0; j < nDim; j++) {
    m_minMax[j] = elements[0].m_minMax[j];
    m_minMax[j + nDim] = elements[0].m_minMax[j + nDim];
  }
  for(MInt i = 0; i < m_noElements; i++) {
    for(MInt j = 0; j < nDim; j++) {
      m_minMax[j + nDim] =
          (m_minMax[j + nDim] < elements[i].m_minMax[j + nDim]) ? elements[i].m_minMax[j + nDim] : m_minMax[j + nDim];
      m_minMax[j] = (m_minMax[j] > elements[i].m_minMax[j]) ? elements[i].m_minMax[j] : m_minMax[j];
    }
  }
  m_log << " Bounding box for segment :" << endl;
  m_log << " max = ";
  for(MInt i = 0; i < nDim; i++) {
    m_log << m_minMax[i + nDim] << " ";
  }
  m_log << endl;
 
  m_log << " min = ";
  for(MInt i = 0; i < nDim; i++) {
    m_log << m_minMax[i] << " ";
  }
  m_log << endl;
 
  if(m_GFieldInitFromSTL) {
    // Calculate bounding box of Moving Boundary segment
    UpdateMBBoundingBox();
    m_log << " Bounding box for Moving Boundary segment:" << endl;
    m_log << " max = ";
    for(MInt i = 0; i < nDim; i++) {
      m_log << m_mbminMax[i + nDim] << " ";
    }
    m_log << endl;
    m_log << " min = ";
    for(MInt i = 0; i < nDim; i++) {
      m_log << m_mbminMax[i] << " ";
    }
    m_log << endl;
  }
}

countSegmentLinesASCII()

void Geometry2D::countSegmentLinesASCII ( const MString &  fileName, MInt *  noElements  )

inlineprotected

Author

Thomas Schilden

Date

27.12.2016

Parameters

[in]	fileName	the name of the file to open
[in]	the	number of entries to be returned (gets nulled in function again)

Definition at line 1012 of file geometry2d.cpp.

                                                                                        {
  TRACE();
 
  *noElements = 0;
 
  m_log << " Counting segment file: " << fileName << endl;
 
  ifstream ifl(fileName);
  if(!ifl) {
    stringstream errorMessage;
    errorMessage << " ERROR in segment::readSegment (counting loop), couldn't find file : " << fileName << "!";
    mTerm(1, AT_, errorMessage.str());
  }
  // If number of elements unknown, count them...
  ifl >> (*noElements);
  ifl.close();
  (*noElements)--; // The ascii files contain the points(!) not the elements
  m_log << " number of vertices = " << *noElements << endl;
}

edgeTriangleIntersection()

MBool Geometry2D::edgeTriangleIntersection ( MFloat *  trianglePoint1, MFloat *  trianglePoint2, MFloat *  trianglePoint3, MFloat *  edgePoint1, MFloat *  edgePoint2  )

inlinevirtual

Reimplemented from Geometry< 2 >.

Definition at line 604 of file geometry2d.cpp.

                                                                      {
  //  TRACE();
 
  MFloat a[2]; // point in plane
  MFloat b[2]; // point in plane
  MFloat c[2]; // start of piercing edge
  MFloat d[2]; // end of piercing edge
  MFloat s1, s2, gamma;
 
 
  // Calculate the intersection between edge and triangle and check if
  // the intersection point lies on the edge.
 
  // pierce plane found
  // TODO labels:GEOM why is a,b,c,d used here?
  for(MInt k = 0; k < nDim; k++) {
    a[k] = edgePoint1[k];
    b[k] = edgePoint2[k];
    c[k] = trianglePoint1[k];
    d[k] = trianglePoint2[k];
  }
 
  gamma = (b[0] - a[0]) * (c[1] - d[1]) - (c[0] - d[0]) * (b[1] - a[1]);
 
  // TODO labels:GEOM,toenhance the small number should be replaced by an eps value defined in maia.h
#ifdef IBM_BLUE_GENE
  if(gamma < 0.00000000001) {
    return false;
  }
#endif
 
  s1 = ((c[0] - d[0]) * (a[1] - c[1]) - (c[1] - d[1]) * (a[0] - c[0])) / gamma;
 
  s2 = ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) / gamma;
 
  if(s2 < 0 || s2 > 1.0 || s1 < 0 || s1 > 1.0) {
    return false;
  } else {
    return true;
  }
}

getIntersectionElements() [1/2]

MInt Geometry2D::getIntersectionElements ( MFloat *  targetRegion, std::vector< MInt > &  nodeList  )

virtual

Author

Rainhill Freitas

Date

unknown

Bug:

labels:GEOM This function cannot be used on IBM_BLUE_GENE. It uses Cramer's rule to solve a 2 by 2 linear equation system to find intersection points. The determinat in the denominator of Cramer's rule can be become 0 which cannot be handled by IBM_BLUE_GENE. Use the other Geometry2D::getIntersectionElements instead.

Reimplemented from Geometry< 2 >.

Definition at line 432 of file geometry2d.cpp.

                                                                                        {
  //  TRACE();
 
  MInt noReallyIntersectingNodes = 0;
  // get all candidates for an intersection
  // Check for intersection...
  bitset<4> points[2];
  bitset<4> faceCodes[4];
  bitset<4> pCode;
 
  // Edges of the targetRegion (using targetPoints)
  MInt rejection;
  MBool piercePointInside;
  MBool triviallyAccepted;
 
  // Each of the following arrays holds one different point for
  // all of the 6 planes. Points are built with the targetRegion
  // i.e. a = targetRegion[pointAInPlane[0]],
  //          targetRegion[pointAInPlane[1]],
  // etc.
  MInt pointAInPlane[4][2] = {{0, 1}, {2, 1}, {0, 1}, {0, 3}};
 
  MInt pointBInPlane[4][2] = {{0, 3}, {2, 3}, {2, 1}, {2, 3}};
 
  MFloat a[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // point in plane
  MFloat b[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // point in plane
  MFloat c[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // start of piercing edge
  MFloat d[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // end of piercing edge
  MFloat s1, s2, gamma;                          // For pierce point calculation
  MFloat pP[2] = {numeric_limits<MFloat>::max(),
                  numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // piercePoint
  faceCodes[0] = IPOW2(0);
  faceCodes[1] = IPOW2(1);
  faceCodes[2] = IPOW2(2);
  faceCodes[3] = IPOW2(3);
  bitset<4> result;
 
  m_adt->retrieveNodes(targetRegion, nodeList);
  const MInt noNodes = nodeList.size();
  //  return noNodes;
  noReallyIntersectingNodes = 0;
 
  for(MInt n = 0; n < noNodes; n++) {
    // create edges (in 3D) AB, AC, BC and point A B C
    // Determine outcode (see Aftosmis, Solution Adaptive Cartesian Grid Methods for Aerodynamic flows ...)
 
    // Loop over all points of an element<2>
    for(MInt p = 0; p < nDim; p++) {
      points[p] = 0;
      // Calculate outcode for point
      for(MInt j = 0; j < nDim; j++) {
        if(elements[nodeList[n]].m_vertices[p][j] < targetRegion[j]) {
          points[p] |= faceCodes[2 * j];
        }
        if(elements[nodeList[n]].m_vertices[p][j] > targetRegion[j + nDim]) {
          points[p] |= faceCodes[2 * j + 1];
        }
      }
      //      m_log << points[p] << endl;
    }
    rejection = 0;
    // check outcode combinations for edges for trivial rejection
    for(MInt i = 0; i < nDim; i++) {
      if((points[0] & points[1]) != 0) {
        rejection++;
        break;
      } else {
        // If one point is inside region the element<2> is trivially accepted
        triviallyAccepted = false;
        for(MInt k = 0; k < nDim; k++) {
          if(points[k] == 0) {
            triviallyAccepted = true;
            rejection = 0;
            break;
          }
        }
        if(triviallyAccepted) {
          break;
        }
        // No trivial rejection, check for rejection of subsegment:
        // For all pierce points!
        // 1. Calculate pierce point:
        //    a - determine plane for pierce point calculation
        //    b - calculate pierce point
        // 2. Check for rejection of new segment
        //    a - calculate new outcode
        //    b - check for containment in pierce planes face
        // 3. If all(!) pierce points are rejected -> reject edge
        // TODO labels:GEOM This algorith might get a problem if a triangle lies completely on
        //    a face.
 
        // 1.a
        result = (points[0] | points[1]);
        piercePointInside = false;
        for(MInt j = 0; j < 2 * nDim; j++) {
          if(result[j] == 1) {
            // pierce plane found
            for(MInt k = 0; k < nDim; k++) {
              a[k] = targetRegion[pointAInPlane[j][k]];
              b[k] = targetRegion[pointBInPlane[j][k]];
              c[k] = elements[nodeList[n]].m_vertices[0][k];
              d[k] = elements[nodeList[n]].m_vertices[1][k];
            }
            gamma = (b[0] - a[0]) * (c[1] - d[1]) - (c[0] - d[0]) * (b[1] - a[1]);
 
            s1 = ((c[0] - d[0]) * (a[1] - c[1]) - (c[1] - d[1]) * (a[0] - c[0])) / gamma;
 
 
            s2 = ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) / gamma;
            // 1. b Pierce point pP in plane j:
            for(MInt k = 0; k < nDim; k++) {
              if(s1 * s1 < s2 * s2)
                pP[k] = c[k] + s2 * (d[k] - c[k]);
              else
                pP[k] = a[k] + s1 * (b[k] - a[k]);
            }
            pCode = 0;
            // 2. a Calculate outcode for pierce point
            for(MInt k = 0; k < nDim; k++) {
              if(pP[k] < targetRegion[k]) {
                pCode |= faceCodes[2 * k];
              }
              if(pP[k] > targetRegion[k + nDim]) {
                pCode |= faceCodes[2 * k + 1];
              }
            }
 
          } else {
            continue;
          }
          // 2. b
          result = faceCodes[j];
          result.flip();
          result = (result & pCode);
          if(result == 0) { // -> is contained
            piercePointInside = true;
            break;
          }
        }
        // reject if all pierce points are off coresponding face
        // else accept
        if(!piercePointInside) {
          rejection++;
        } else {
          rejection = 0;
          break;
        }
      }
    }
    // If not all edges are rejected a cutting element<2> has been found
    //    m_log << rejection << endl;
    if(!rejection) {
      nodeList[noReallyIntersectingNodes] = nodeList[n];
      noReallyIntersectingNodes++;
      continue;
    }
 
    // write nodes in reallyIntersectingNodes
  }
  nodeList.resize(noReallyIntersectingNodes);
 
  return noReallyIntersectingNodes;
}

getIntersectionElements() [2/2]

MInt Geometry2D::getIntersectionElements ( MFloat *  targetRegion, std::vector< MInt > &  nodeList, MFloat  cellHalfLength, const MFloat *const  cell_coords  )

virtual

Author

Andreas Lintermann

Date

09.11.2009

The separating axis theorem (SAT) states, that two convex polyhedra, that are to be tested against intersection are disjoint, if they can be separated along either an axis parallel to a normal of a face of the first or the second polyhedra, or along an axis formed from the cross product of an edge from them.

Algorithm:

• 1) Test projection onto \(x,y\)-axis Project the triangle and the gridcell onto both axes and check overlapping

• 2) Test if we have parallel lines. If this is the case we are done

• 3) Otherwise, rotate everything so that the axis to be tested becomes axis-aligned. Then test again agains projection overlapping.

This function replaces the other Geometry2D::getIntersectionElements, which cannot be used on IBM_BLUE_GENE.

Reimplemented from Geometry< 2 >.

Definition at line 310 of file geometry2d.cpp.

                                                                {
  //  TRACE();
 
 
  MInt noReallyIntersectingNodes = 0;
  m_adt->retrieveNodes(targetRegion, nodeList);
  const MInt noNodes = nodeList.size();
 
  MFloat cos_alpha, sin_alpha, length_diff;
  MFloat diff_vec[nDim];
 
  for(MInt i = 0; i < nDim; i++) {
    diff_vec[i] = numeric_limits<MFloat>::max(); // gcc 4.8.2 maybe uninitialized
  }
 
  MFloat rot_point;
  MFloat rot_corners[nDim];
 
  MBool overlap[2];
  MBool parallel;
  MInt parallel_dim;
 
  for(MInt n = 0; n < noNodes; n++) {
    length_diff = 0.0;
    parallel = false;
    overlap[0] = false;
    overlap[1] = false;
 
    parallel_dim = 0;
 
    // coordinate-axis tests
    for(MInt i = 0; i < nDim; i++) // run over dimensions (x and y)
    {
      // is first point in projection of quad?
      if(elements[nodeList[n]].m_vertices[0][i] >= targetRegion[i]
         && elements[nodeList[n]].m_vertices[0][i] <= targetRegion[i + nDim]) {
        overlap[i] = true;
        continue;
      }
 
      // is second point in projection of quad?
      if(elements[nodeList[n]].m_vertices[1][i] >= targetRegion[i]
         && elements[nodeList[n]].m_vertices[1][i] <= targetRegion[i + nDim]) {
        overlap[i] = true;
        continue;
      }
 
      // do both points clasp the projection of quad?
      if(elements[nodeList[n]].m_vertices[0][i] <= targetRegion[i]
         && elements[nodeList[n]].m_vertices[1][i] >= targetRegion[i + nDim]) {
        overlap[i] = true;
        continue;
      }
      if(elements[nodeList[n]].m_vertices[1][i] <= targetRegion[i]
         && elements[nodeList[n]].m_vertices[0][i] >= targetRegion[i + nDim]) {
        overlap[i] = true;
        continue;
      }
      if(!overlap[i]) break;
    }
 
    if(!(overlap[0] && overlap[1])) continue;
 
    // is the line parallel to axis? Then we are done if we don't
    for(parallel_dim = 0; parallel_dim < nDim; parallel_dim++) {
      parallel = parallel
                 || (approx(elements[nodeList[n]].m_vertices[0][parallel_dim],
                            elements[nodeList[n]].m_vertices[1][parallel_dim], MFloatEps));
      if(parallel) break;
    }
 
    if(parallel) {
      nodeList[noReallyIntersectingNodes] = nodeList[n];
      noReallyIntersectingNodes++;
    } else {
      // For the sake of simplicity rotate everything so that line becomes parallel to x-axis. Check overlapping in
      // y-direction.
      for(MInt i = 0; i < nDim; i++) {
        diff_vec[i] = elements[nodeList[n]].m_vertices[1][i] - elements[nodeList[n]].m_vertices[0][i];
        length_diff += diff_vec[i] * diff_vec[i];
      }
      length_diff = sqrt(length_diff);
 
      cos_alpha = fabs(diff_vec[0]) / length_diff;
      sin_alpha = fabs(diff_vec[1]) / length_diff;
 
      if((diff_vec[0] * diff_vec[1]) < 0) // negative slope
      {
        rot_point =
            sin_alpha * elements[nodeList[n]].m_vertices[0][0] + cos_alpha * elements[nodeList[n]].m_vertices[0][1];
        rot_corners[0] = sin_alpha * (targetRegion[2]) + cos_alpha * (targetRegion[3]);
        rot_corners[1] = sin_alpha * (targetRegion[0]) + cos_alpha * (targetRegion[1]);
      } else // positive slope
      {
        rot_point = -1 * sin_alpha * elements[nodeList[n]].m_vertices[0][0]
                    + cos_alpha * elements[nodeList[n]].m_vertices[0][1];
        rot_corners[0] = -1 * sin_alpha * (targetRegion[0]) + cos_alpha * (targetRegion[3]);
        rot_corners[1] = -1 * sin_alpha * (targetRegion[2]) + cos_alpha * (targetRegion[1]);
      }
 
      // final cut test
      if(rot_corners[0] >= rot_point && rot_corners[1] <= rot_point) {
        nodeList[noReallyIntersectingNodes] = nodeList[n];
        noReallyIntersectingNodes++;
      }
    }
  }
  nodeList.resize(noReallyIntersectingNodes);
 
  return noReallyIntersectingNodes;
}

getIntersectionMBElements()

MInt Geometry2D::getIntersectionMBElements ( MFloat *  targetRegion, std::vector< MInt > &  nodeList  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 688 of file geometry2d.cpp.

                                                                                          {
  //   TRACE();
 
  MInt noReallyIntersectingNodes = 0;
  // get all candidates for an intersection
  // Check for intersection...
  bitset<4> points[2];
  bitset<4> faceCodes[4];
  bitset<4> pCode;
 
  // Edges of the targetRegion (using targetPoints)
  MInt rejection;
  MBool piercePointInside;
  MBool triviallyAccepted;
 
  // Each of the following arrays holds one different point for
  // all of the 6 planes. Points are built with the targetRegion
  // i.e. a = targetRegion[pointAInPlane[0]],
  //          targetRegion[pointAInPlane[1]],
  // etc.
  MInt pointAInPlane[4][2] = {{0, 1}, {2, 1}, {0, 1}, {0, 3}};
 
  MInt pointBInPlane[4][2] = {{0, 3}, {2, 3}, {2, 1}, {2, 3}};
 
  MFloat a[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // point in plane
  MFloat b[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // point in plane
  MFloat c[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // start of piercing edge
  MFloat d[2] = {numeric_limits<MFloat>::max(),
                 numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // end of piercing edge
  MFloat s1, s2, gamma;                          // For pierce point calculation
  MFloat pP[2] = {numeric_limits<MFloat>::max(),
                  numeric_limits<MFloat>::max()}; // gcc 4.8.2 maybe uninitialized // piercePoint
  faceCodes[0] = IPOW2(0);
  faceCodes[1] = IPOW2(1);
  faceCodes[2] = IPOW2(2);
  faceCodes[3] = IPOW2(3);
  bitset<4> result;
 
  m_adt->retrieveNodesMBElements(targetRegion, nodeList);
  const MInt noNodes = nodeList.size();
  //  return noNodes;
  noReallyIntersectingNodes = 0;
 
  for(MInt n = 0; n < noNodes; n++) {
    // create edges (in 3D) AB, AC, BC and point A B C
    // Determine outcode (see Aftosmis, Solution Adaptive Cartesian Grid Methods for Aerodynamic flows ...)
 
    // Loop over all points of an element<2>
    for(MInt p = 0; p < nDim; p++) {
      points[p] = 0;
      // Calculate outcode for point
      for(MInt j = 0; j < nDim; j++) {
        if(mbelements[nodeList[n]].m_vertices[p][j] < targetRegion[j]) {
          points[p] |= faceCodes[2 * j];
        }
        if(mbelements[nodeList[n]].m_vertices[p][j] > targetRegion[j + nDim]) {
          points[p] |= faceCodes[2 * j + 1];
        }
      }
      //      m_log << points[p] << endl;
    }
    rejection = 0;
    // check outcode combinations for edges for trivial rejection
    for(MInt i = 0; i < nDim; i++) {
      if((points[0] & points[1]) != 0) {
        rejection++;
        break;
      } else {
        // If one point is inside region the element<2> is trivially accepted
        triviallyAccepted = false;
        for(MInt k = 0; k < nDim; k++) {
          if(points[k] == 0) {
            triviallyAccepted = true;
            rejection = 0;
            break;
          }
        }
        if(triviallyAccepted) {
          break;
        }
        // No trivial rejection, check for rejection of subsegment:
        // For all pierce points!
        // 1. Calculate pierce point:
        //    a - determine plane for pierce point calculation
        //    b - calculate pierce point
        // 2. Check for rejection of new segment
        //    a - calculate new outcode
        //    b - check for containment in pierce planes face
        // 3. If all(!) pierce points are rejected -> reject edge
        // TODO labels:GEOM This algorith might get a problem if a triangle lies completely on
        //    a face.
 
        // 1.a
        result = (points[0] | points[1]);
        piercePointInside = false;
        for(MInt j = 0; j < 2 * nDim; j++) {
          if(result[j] == 1) {
            // pierce plane found
            for(MInt k = 0; k < nDim; k++) {
              a[k] = targetRegion[pointAInPlane[j][k]];
              b[k] = targetRegion[pointBInPlane[j][k]];
              c[k] = mbelements[nodeList[n]].m_vertices[0][k];
              d[k] = mbelements[nodeList[n]].m_vertices[1][k];
            }
            gamma = (b[0] - a[0]) * (c[1] - d[1]) - (c[0] - d[0]) * (b[1] - a[1]);
 
            s1 = ((c[0] - d[0]) * (a[1] - c[1]) - (c[1] - d[1]) * (a[0] - c[0])) / gamma;
 
 
            s2 = ((b[0] - a[0]) * (c[1] - a[1]) - (c[0] - a[0]) * (b[1] - a[1])) / gamma;
            // 1. b Pierce point pP in plane j:
            for(MInt k = 0; k < nDim; k++) {
              if(s1 * s1 < s2 * s2)
                pP[k] = c[k] + s2 * (d[k] - c[k]);
              else
                pP[k] = a[k] + s1 * (b[k] - a[k]);
            }
            pCode = 0;
            // 2. a Calculate outcode for pierce point
            for(MInt k = 0; k < nDim; k++) {
              if(pP[k] < targetRegion[k]) {
                pCode |= faceCodes[2 * k];
              }
              if(pP[k] > targetRegion[k + nDim]) {
                pCode |= faceCodes[2 * k + 1];
              }
            }
 
          } else {
            continue;
          }
          // 2. b
          result = faceCodes[j];
          result.flip();
          result = (result & pCode);
          if(result == 0) { // -> is contained
            piercePointInside = true;
            break;
          }
        }
        // reject if all pierce points are off coresponding face
        // else accept
        if(!piercePointInside) {
          rejection++;
        } else {
          rejection = 0;
          break;
        }
      }
    }
    // If not all edges are rejected a cutting element<2> has been found
    //    m_log << rejection << endl;
    if(!rejection) {
      nodeList[noReallyIntersectingNodes] = nodeList[n];
      noReallyIntersectingNodes++;
      continue;
    }
 
    // write nodes in reallyIntersectingNodes
  }
  nodeList.resize(noReallyIntersectingNodes);
 
  return noReallyIntersectingNodes;
}

getLineIntersectionElements()

MInt Geometry2D::getLineIntersectionElements ( MFloat *  targetRegion, std::vector< MInt > &  nodeList  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 654 of file geometry2d.cpp.

                                                                                            {
  //  TRACE();
 
  MInt noReallyIntersectingNodes = 0;
 
  m_adt->retrieveNodes(targetRegion, nodeList);
  const MInt noNodes = nodeList.size();
 
  MFloat e1[2], e2[2];
  for(MInt i = 0; i < nDim; i++) {
    e1[i] = targetRegion[i];
    e2[i] = targetRegion[i + nDim];
  }
 
  for(MInt n = 0; n < noNodes; n++) {
    if(edgeTriangleIntersection(elements[nodeList[n]].m_vertices[0], elements[nodeList[n]].m_vertices[1], 0, e1, e2)) {
      nodeList[noReallyIntersectingNodes] = nodeList[n];
      noReallyIntersectingNodes++;
    }
  }
  nodeList.resize(noReallyIntersectingNodes);
 
  return noReallyIntersectingNodes;
}

getLineIntersectionElementsOld1()

MInt Geometry2D::getLineIntersectionElementsOld1 ( MFloat *  targetRegion, std::vector< MInt > &  nodeList  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 679 of file geometry2d.cpp.

                                                                                                {
  return getLineIntersectionElements(targetRegion, nodeList);
}

getLineIntersectionElementsOld2()

MInt Geometry2D::getLineIntersectionElementsOld2 ( MFloat *  targetRegion, MInt *  spaceDirection, std::vector< MInt > &  nodeList  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 683 of file geometry2d.cpp.

                                                                            {
  return getLineIntersectionElements(targetRegion, nodeList);
}

getLineIntersectionMBElements()

MInt Geometry2D::getLineIntersectionMBElements ( MFloat *  targetRegion, std::vector< MInt > &  nodeList  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 856 of file geometry2d.cpp.

                                                                                              {
  //   TRACE();
 
  MInt noReallyIntersectingNodes = 0;
 
  m_adt->retrieveNodesMBElements(targetRegion, nodeList);
  const MInt noNodes = nodeList.size();
 
  MFloat e1[2], e2[2];
  for(MInt i = 0; i < nDim; i++) {
    e1[i] = targetRegion[i];
    e2[i] = targetRegion[i + nDim];
  }
 
  for(MInt n = 0; n < noNodes; n++) {
    if(edgeTriangleIntersection(mbelements[nodeList[n]].m_vertices[0], mbelements[nodeList[n]].m_vertices[1], 0, e1,
                                e2)) {
      nodeList[noReallyIntersectingNodes] = nodeList[n];
      noReallyIntersectingNodes++;
    }
  }
  nodeList.resize(noReallyIntersectingNodes);
 
  return noReallyIntersectingNodes;
}

getSphereIntersectionMBElements()

MInt Geometry2D::getSphereIntersectionMBElements ( MFloat *  P, MFloat  radius, std::vector< MInt > &  nodeList  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 883 of file geometry2d.cpp.

                                                                                                    {
  // TRACE();
 
  MInt noReallyIntersectingNodes = 0;
 
  // compute minimum target region that contains the sphere:
  MFloat target[4] = {0, 0, 0, 0};
  MFloat enlargeFactor = 1.05;
  for(MInt i = 0; i < nDim; i++) {
    target[i] = P[i] - radius * enlargeFactor;
    target[i + nDim] = P[i] + radius * enlargeFactor;
  }
 
  // fetch all possibly intersecting triangles in this region:
  m_adt->retrieveNodesMBElements(target, nodeList);
  const MInt noNodes = nodeList.size();
 
  MFloat A[2], B[2], AB[2], Q1[2];
 
  // loop over all elements
  for(MInt n = 0; n < noNodes; n++) {
    // store the vertices of the current element<2> in A, B, C:
    for(MInt i = 0; i < nDim; i++) {
      A[i] = mbelements[nodeList[n]].m_vertices[0][i];
      B[i] = mbelements[nodeList[n]].m_vertices[1][i];
    }
 
    // compute separation algorithm from http://realtimecollisiondetection.net/blog/?p=103:
    for(MInt i = 0; i < nDim; i++) {
      A[i] = A[i] - P[i];
      B[i] = B[i] - P[i];
      AB[i] = B[i] - A[i];
    }
    MFloat rr = radius * radius;
    MFloat aa = F0;
    MFloat ab = F0;
    MFloat bb = F0;
    MFloat e1 = F0;
    for(MInt i = 0; i < nDim; i++) {
      aa += A[i] * A[i];
      ab += A[i] * B[i];
      bb += B[i] * B[i];
      e1 += AB[i] * AB[i];
    }
    MBool sep2 = (aa > rr) && (ab > aa);
    MBool sep3 = (bb > rr) && (ab > bb);
    MFloat d1 = ab - aa;
    MFloat qq1 = F0;
    for(MInt i = 0; i < nDim; i++) {
      Q1[i] = A[i] * e1 - d1 * AB[i];
      qq1 += Q1[i] * Q1[i];
    }
    MBool sep5 = (qq1 > rr * e1 * e1);
 
    MBool separated = sep2 || sep3 || sep5;
 
    if(!separated) {
      nodeList[noReallyIntersectingNodes] = nodeList[n];
      noReallyIntersectingNodes++;
    }
  }
  nodeList.resize(noReallyIntersectingNodes);
 
  return noReallyIntersectingNodes;
}

MoveAllMBElementVertex()

void Geometry2D::MoveAllMBElementVertex ( MFloat *  dx )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 950 of file geometry2d.cpp.

                                                  {
  TRACE();
 
  for(MInt e = 0; e < m_noMBElements; e++) {
    for(MInt j = 0; j < nDim; j++) {
      for(MInt i = 0; i < nDim; i++) {
        mbelements[e].m_vertices[j][i] += dx[i];
      }
    }
  }
}

MoveMBElementVertex()

void Geometry2D::MoveMBElementVertex ( MInt  e, MInt  v, MFloat *  dx  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 962 of file geometry2d.cpp.

                                                               {
  TRACE();
 
  for(MInt i = 0; i < nDim; i++) {
    mbelements[e].m_vertices[v][i] += dx[i];
  }
}

readSegmentLinesASCII()

void Geometry2D::readSegmentLinesASCII ( MString  fileName, Collector< element< 2 > > *  elemCollector, MInt  bndCndId, MInt  segmentId, MInt *  offset  )

inlineprotected

Author

Thomas Schilden

Date

27.12.2016

Parameters

[in]	fileName	the name of the file to open
[in]	etc..

Definition at line 1041 of file geometry2d.cpp.

                                                                            {
  TRACE();
  m_log << " Reading segment file: " << fileName << endl;
  m_log << " BC : " << bndCndId << endl;
 
 
  ifstream ifl(fileName);
  if(!ifl) {
    stringstream errorMessage;
    errorMessage << " ERROR in segment::readSegment (reading loop), couldn't find file : " << fileName << "!";
    mTerm(1, AT_, errorMessage.str());
  }
 
  element<2>* allelements = elemCollector->a;
 
  // Read no of coordinates
  MInt tmp;
  ifl >> tmp;
  elemCollector->append();
  for(MInt j = 0; j < 2; j++) {
    ifl >> allelements[*offset].m_vertices[0][j];
    allelements[*offset].m_normal[j] = F0;
  }
  MInt fileElementCounter = 0;
  while(!ifl.eof()) {
    // Read coordinates
    for(MInt j = 0; j < 2; j++) {
      ifl >> allelements[*offset].m_vertices[1][j];
      allelements[*offset].m_normal[j] = F0;
    }
    allelements[*offset].boundingBox();
    allelements[*offset].m_bndCndId = bndCndId;
    allelements[*offset].m_segmentId = segmentId;
    if(m_GFieldInitFromSTL) {
      for(MInt id = 0; id < m_noLevelSetIntfBndIds; id++) {
        if(bndCndId == m_levelSetIntfBndIds[id]) {
          m_noMBElements++;
        }
      }
    }
 
    if(m_noElements > 1000 && !((*offset) % (m_noElements / 10))) {
      m_log << (MInt)(((*offset) - 1) / ((MFloat)m_noElements / 100)) + 1 << " % read." << endl;
    }
    fileElementCounter++;
    (*offset)++;
    // Check if already all points read
    if(fileElementCounter < tmp - 1) { // set start point of next line
      elemCollector->append();
      for(MInt j = 0; j < 2; j++) {
        allelements[*offset].m_vertices[0][j] = allelements[(*offset) - 1].m_vertices[1][j];
        allelements[*offset].m_normal[j] = F0;
      }
    } else { // If all points read leave loop
      break;
    }
  }
  ifl.close();
}

readSegments()

void Geometry2D::readSegments ( )

protectedvirtual

Reimplemented from Geometry< 2 >.

Definition at line 141 of file geometry2d.cpp.

                              {
  TRACE();
 
  MInt counter = 0;
  MInt segmentId = 0;
  MString fileName;
  // This loop only counts all elements of all segments
  for(m_bodyIt = m_bodyMap.begin(); m_bodyIt != m_bodyMap.end(); m_bodyIt++) {
    // Do not create default body!
    if(m_bodyIt->second->name != "default") {
      for(MInt i = 0; i < m_bodyIt->second->noSegments; i++) {
        fileName = *geometryContext().getProperty("filename", segmentId)->asString();
        if(geometryContext().noPropertySegments("filename") == 1
           && (m_bodyIt->second->noSegments > 1 || m_bodyMap.size() > 2)) {
          fileName += "." + to_string(segmentId);
        }
        countSegmentLinesASCII(fileName, &counter);
        m_noElements += counter;
        segmentId++;
      }
    }
  }
 
  auto* m_allelements = new Collector<element<nDim>>(m_noElements, nDim, 0);
  element<nDim>* allelements = m_allelements->a;
  vector<MInt> allBCs;
 
  segmentId = 0;
  counter = 0;
  // This loop reads all elements from all segments
  for(m_bodyIt = m_bodyMap.begin(); m_bodyIt != m_bodyMap.end(); m_bodyIt++) {
    // Do not create default body!
    if(m_bodyIt->second->name != "default") {
      for(MInt i = 0; i < m_bodyIt->second->noSegments; i++) {
        m_log << " Reading Ascii coordinates from file " << endl;
        fileName = *geometryContext().getProperty("filename", segmentId)->asString();
        if(geometryContext().noPropertySegments("filename") == 1
           && (m_bodyIt->second->noSegments > 1 || m_bodyMap.size() > 2)) {
          fileName += "." + to_string(segmentId);
        }
        MInt bndCndId = *geometryContext().getProperty("BC", segmentId)->asInt();
        allBCs.push_back(bndCndId);
        readSegmentLinesASCII(fileName, m_allelements, bndCndId, segmentId, &counter);
        segmentId++;
      }
    }
  }
  m_log << m_noElements << " line elements read. " << endl;
 
  if(m_GFieldInitFromSTL) {
    m_mbelements = new Collector<element<nDim>>(m_noMBElements, nDim, 0);
    m_noElements -= m_noMBElements;
    mbelements = m_mbelements->a;
  }
 
  m_elements = new Collector<element<nDim>>(m_noElements, nDim, 0);
  elements = m_elements->a;
 
  MInt mbelem_counter = 0, elem_counter = 0;
  MBool mbElem = false;
 
  for(MInt allelem_counter = 0; allelem_counter < m_allelements->size(); allelem_counter++) {
    if(m_GFieldInitFromSTL && m_noLevelSetIntfBndIds > 0) {
      for(MInt id = 0; id < m_noLevelSetIntfBndIds; id++) {
        if(allelements[allelem_counter].m_bndCndId == m_levelSetIntfBndIds[id]) {
          mbElem = true;
          break;
        }
      }
      if(mbElem) {
        // LevelSet Moving boundary elements
        m_mbelements->append();
        for(MInt j = 0; j < 2; j++) {
          mbelements[mbelem_counter].m_normal[j] = allelements[allelem_counter].m_normal[j];
          for(MInt i = 0; i < 2; i++) {
            mbelements[mbelem_counter].m_vertices[j][i] = allelements[allelem_counter].m_vertices[j][i];
          }
        }
        mbelements[mbelem_counter].boundingBox();
        mbelements[mbelem_counter].m_bndCndId = allelements[allelem_counter].m_bndCndId;
        mbelements[mbelem_counter].m_segmentId = allelements[allelem_counter].m_segmentId;
        mbelem_counter++;
        mbElem = false;
      } else {
        // Non movable elements
        m_elements->append();
        for(MInt j = 0; j < 2; j++) {
          elements[elem_counter].m_normal[j] = allelements[allelem_counter].m_normal[j];
          for(MInt i = 0; i < 2; i++) {
            elements[elem_counter].m_vertices[j][i] = allelements[allelem_counter].m_vertices[j][i];
          }
        }
        elements[elem_counter].boundingBox();
        elements[elem_counter].m_bndCndId = allelements[allelem_counter].m_bndCndId;
        elements[elem_counter].m_segmentId = allelements[allelem_counter].m_segmentId;
        elem_counter++;
      }
    } else {
      if(m_GFieldInitFromSTL && (m_levelSetIntfBndId == allelements[allelem_counter].m_bndCndId)) {
        // LevelSet Moving boundary elements
        m_mbelements->append();
        for(MInt i = 0; i < 2; i++) {
          mbelements[mbelem_counter].m_normal[i] = allelements[allelem_counter].m_normal[i];
        }
        for(MInt j = 0; j < 2; j++) {
          for(MInt i = 0; i < 2; i++) {
            mbelements[mbelem_counter].m_vertices[j][i] = allelements[allelem_counter].m_vertices[j][i];
          }
        }
        mbelements[mbelem_counter].boundingBox();
        mbelements[mbelem_counter].m_bndCndId = allelements[allelem_counter].m_bndCndId;
        mbelements[mbelem_counter].m_segmentId = allelements[allelem_counter].m_segmentId;
        mbelem_counter++;
      } else {
        // Non movable elements
        m_elements->append();
        for(MInt i = 0; i < 2; i++) {
          elements[elem_counter].m_normal[i] = allelements[allelem_counter].m_normal[i];
        }
        for(MInt j = 0; j < 2; j++) {
          for(MInt i = 0; i < 2; i++) {
            elements[elem_counter].m_vertices[j][i] = allelements[allelem_counter].m_vertices[j][i];
          }
        }
        elements[elem_counter].boundingBox();
        elements[elem_counter].m_bndCndId = allelements[allelem_counter].m_bndCndId;
        elements[elem_counter].m_segmentId = allelements[allelem_counter].m_segmentId;
        elem_counter++;
      }
    }
  }
 
  delete m_allelements;
 
  calculateBoundingBox();
 
  m_noAllBCs = allBCs.size();
  mAlloc(m_allBCs, m_noAllBCs, "m_allBCs", 0, AT_);
  for(MInt i = 0; i < m_noAllBCs; i++)
    m_allBCs[i] = allBCs[i];
}

ReplaceMBElementVertex()

void Geometry2D::ReplaceMBElementVertex ( MInt  e, MInt  v, MFloat *  np  )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 970 of file geometry2d.cpp.

                                                                  {
  TRACE();
 
  for(MInt i = 0; i < nDim; i++)
    mbelements[e].m_vertices[v][i] = np[i];
}

UpdateADT()

void Geometry2D::UpdateADT ( )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 997 of file geometry2d.cpp.

                           {
  TRACE();
 
  m_adt->buildTreeMB();
}

UpdateMBBoundingBox()

void Geometry2D::UpdateMBBoundingBox ( )

virtual

Reimplemented from Geometry< 2 >.

Definition at line 977 of file geometry2d.cpp.

                                     {
  TRACE();
 
  for(MInt e = 0; e < m_mbelements->size(); e++) {
    mbelements[e].boundingBox();
  }
  for(MInt j = 0; j < nDim; j++) {
    m_mbminMax[j] = mbelements[0].m_minMax[j];
    m_mbminMax[j + nDim] = mbelements[0].m_minMax[j + nDim];
  }
  for(MInt i = 0; i < m_mbelements->size(); i++) {
    for(MInt j = 0; j < nDim; j++) {
      m_mbminMax[j + nDim] = (m_mbminMax[j + nDim] < mbelements[i].m_minMax[j + nDim])
                                 ? mbelements[i].m_minMax[j + nDim]
                                 : m_mbminMax[j + nDim];
      m_mbminMax[j] = (m_mbminMax[j] > mbelements[i].m_minMax[j]) ? mbelements[i].m_minMax[j] : m_mbminMax[j];
    }
  }
}

UpdateMBNormalVector()

virtual void Geometry2D::UpdateMBNormalVector ( MInt  )

inlinevirtual

Reimplemented from Geometry< 2 >.

Definition at line 36 of file geometry2d.h.

{};

Member Data Documentation

m_GFieldInitFromSTL

MBool Geometry2D::m_GFieldInitFromSTL

protected

Definition at line 47 of file geometry2d.h.

m_levelSetIntfBndId

MInt Geometry2D::m_levelSetIntfBndId {}

protected

Definition at line 48 of file geometry2d.h.

m_levelSetIntfBndIds

MInt* Geometry2D::m_levelSetIntfBndIds {}

protected

Definition at line 49 of file geometry2d.h.

m_noLevelSetIntfBndIds

MInt Geometry2D::m_noLevelSetIntfBndIds {}

protected

Definition at line 50 of file geometry2d.h.

nDim

constexpr const MInt Geometry2D::nDim = 2

staticconstexprprotected

Definition at line 52 of file geometry2d.h.

---

The documentation for this class was generated from the following files:

• /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/GEOM/geometry2d.h
• /home/gitlab-runner/scratch/builds/oxpnswJ6/1/aia/m-AIA/m-AIA/src/GEOM/geometry2d.cpp