[cig-commits] r17219 - in short/3D/PyLith/trunk: libsrc/topology unittests/libtests/topology unittests/libtests/topology/data
brad at geodynamics.org
brad at geodynamics.org
Fri Sep 24 17:58:45 PDT 2010
Author: brad
Date: 2010-09-24 17:58:45 -0700 (Fri, 24 Sep 2010)
New Revision: 17219
Modified:
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc
short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh
short/3D/PyLith/trunk/unittests/libtests/topology/data/MeshDataCohesiveTri3Level2Fault1.cc
Log:
First draft of reimplementation of global refinement. Still need to use template for CellRefiner and refactor updating labels.
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.cc 2010-09-25 00:58:45 UTC (rev 17219)
@@ -110,8 +110,9 @@
break;
} // TRIANGLE
case LINE_COHESIVE_LAGRANGE: {
- const int coneVertexOffset = orderNewMesh.verticesCensored().min() - orderOldMesh.verticesCensored().min();
- _newCells_LINE_COHESIVE_LAGRANGE(cells, numCells, cone, coneSize, coneVertexOffset);
+ const int coneVertexOffsetNormal = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
+ const int coneVertexOffsetCensored = orderNewMesh.verticesCensored().min() - orderOldMesh.verticesCensored().min();
+ _newCells_LINE_COHESIVE_LAGRANGE(cells, numCells, cone, coneSize, coneVertexOffsetNormal, coneVertexOffsetCensored);
break;
} // LINE_COHESIVE_LAGRANGE
default:
@@ -154,6 +155,53 @@
} // setCoordsNewVertices
// ----------------------------------------------------------------------
+// Add space for new vertices in group.
+void
+ALE::CellRefinerTri3::groupAddNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup)
+{ // groupAddNewVertices
+ assert(!newGroup.isNull());
+ assert(!oldGroup.isNull());
+
+ const edge_map_type::const_iterator edgesEnd = _edgeToVertex.end();
+ for (edge_map_type::const_iterator e_iter = _edgeToVertex.begin(); e_iter != edgesEnd; ++e_iter) {
+ const point_type newVertex = e_iter->second;
+ const point_type edgeVertexA = e_iter->first.first;
+ const point_type edgeVertexB = e_iter->first.second;
+
+ if (oldGroup->getFiberDimension(edgeVertexA) && oldGroup->getFiberDimension(edgeVertexB)) {
+ if (oldGroup->restrictPoint(edgeVertexA)[0] == oldGroup->restrictPoint(edgeVertexB)[0]) {
+ newGroup->setFiberDimension(newVertex, 1);
+ } // if
+ } // if
+ } // for
+} // groupAddNewVertices
+
+// ----------------------------------------------------------------------
+// Set new vertices in group.
+void
+ALE::CellRefinerTri3::groupSetNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup)
+{ // groupSetNewVertices
+ assert(!newGroup.isNull());
+ assert(!oldGroup.isNull());
+
+ const edge_map_type::const_iterator edgesEnd = _edgeToVertex.end();
+ for (edge_map_type::const_iterator e_iter = _edgeToVertex.begin(); e_iter != edgesEnd; ++e_iter) {
+ const point_type newVertex = e_iter->second;
+ const point_type edgeVertexA = e_iter->first.first;
+ const point_type edgeVertexB = e_iter->first.second;
+
+ if (oldGroup->getFiberDimension(edgeVertexA) && oldGroup->getFiberDimension(edgeVertexB)) {
+ if (oldGroup->restrictPoint(edgeVertexA)[0] == oldGroup->restrictPoint(edgeVertexB)[0]) {
+ newGroup->updatePoint(newVertex, oldGroup->restrictPoint(edgeVertexA));
+ std::cout << "Adding new vertex: " << newVertex << " based on old vertices " << edgeVertexA << " and " << edgeVertexB << std::endl;
+ } // if
+ } // if
+ } // for
+} // groupSetNewVertices
+
+// ----------------------------------------------------------------------
// Get cell type.
ALE::CellRefinerTri3::CellEnum
ALE::CellRefinerTri3::_cellType(const point_type cell)
@@ -276,7 +324,8 @@
int *numCells,
const point_type cone[],
const int coneSize,
- const int coneVertexOffset)
+ const int coneVertexOffsetNormal,
+ const int coneVertexOffsetCensored)
{ // _newCells_LINE_COHESIVE_LAGRANGE
const int coneSizeLine6 = 6;
const int numEdgesLine6 = 3;
@@ -298,23 +347,24 @@
} // for
// new cell 0
- lineCells[0*6+0] = cone[0] + coneVertexOffset;
+ lineCells[0*6+0] = cone[0] + coneVertexOffsetNormal;
lineCells[0*6+1] = newVertices[0];
- lineCells[0*6+2] = cone[2] + coneVertexOffset;
+ lineCells[0*6+2] = cone[2] + coneVertexOffsetNormal;
lineCells[0*6+3] = newVertices[1];
- lineCells[0*6+4] = cone[4] + coneVertexOffset;
+ lineCells[0*6+4] = cone[4] + coneVertexOffsetCensored;
lineCells[0*6+5] = newVertices[2];
// new cell 1
lineCells[1*6+0] = newVertices[0];
- lineCells[1*6+1] = cone[1] + coneVertexOffset;
+ lineCells[1*6+1] = cone[1] + coneVertexOffsetNormal;
lineCells[1*6+2] = newVertices[1];
- lineCells[1*6+3] = cone[3] + coneVertexOffset;
+ lineCells[1*6+3] = cone[3] + coneVertexOffsetNormal;
lineCells[1*6+4] = newVertices[2];
- lineCells[1*6+5] = cone[5] + coneVertexOffset;
+ lineCells[1*6+5] = cone[5] + coneVertexOffsetCensored;
*numCells = 2;
*cells = lineCells;
} // _newCells_LINE_COHESIVE_LAGRANGE
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerTri3.hh 2010-09-25 00:58:45 UTC (rev 17219)
@@ -94,6 +94,22 @@
void setCoordsNewVertices(const ALE::Obj<mesh_type::real_section_type>& newCoordsSection,
const ALE::Obj<mesh_type::real_section_type>& oldCoordsSection);
+ /** Add space for new vertices in group.
+ *
+ * @param newGroup Group in refine mesh.
+ * @param oldGroup Group in original mesh.
+ */
+ void groupAddNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup);
+
+ /** Set new vertices in group.
+ *
+ * @param newGroup Group in refine mesh.
+ * @param oldGroup Group in original mesh.
+ */
+ void groupSetNewVertices(const ALE::Obj<mesh_type::int_section_type>& newGroup,
+ const ALE::Obj<mesh_type::int_section_type>& oldGroup);
+
// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
private :
@@ -176,13 +192,15 @@
* @param numCells Number of refined cells.
* @param cone Vertices in cell (original mesh).
* @param coneSize Number of vertices in cell.
- * @param coneVertexOffset Offset for cone vertices.
+ * @param coneVertexOffsetNormal Offset for normal cone vertices.
+ * @param coneVertexOffset Offset for censored cone vertices.
*/
void _newCells_LINE_COHESIVE_LAGRANGE(const point_type** cells,
int *numCells,
const point_type cone[],
const int coneSize,
- const int coneVertexOffset);
+ const int coneVertexOffsetNormal,
+ const int coneVertexOffsetCensored);
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshOrder.cc 2010-09-25 00:58:45 UTC (rev 17219)
@@ -39,19 +39,48 @@
void
ALE::MeshOrder::initialize(const ALE::Obj<mesh_type>& mesh)
{ // initialize
+ assert(!mesh.isNull());
+
const ALE::Obj<mesh_type::label_sequence>& cells = mesh->heightStratum(0);
assert(!cells.isNull());
- const int numCells = cells->size();
- _cellsNormal = ALE::Interval<point_type>(0, numCells);
-
+
const Obj<mesh_type::label_sequence>& vertices = mesh->depthStratum(0);
assert(!vertices.isNull());
- const int numVertices = vertices->size();
- _verticesNormal = ALE::Interval<point_type>(numCells, numCells+numVertices);
- const int numEntities = numCells + numVertices;
- _cellsCensored = ALE::Interval<point_type>(numEntities, numEntities);
- _verticesCensored = ALE::Interval<point_type>(numEntities, numEntities);
+ if (mesh->hasLabel("censored depth")) {
+ // Count number of cells in censored depth (normal cells).
+ const Obj<mesh_type::label_sequence>& cellsNormal = mesh->getLabelStratum("censored depth", mesh->depth());
+ assert(!cellsNormal.isNull());
+ const mesh_type::label_sequence::iterator cellsNormalEnd = cellsNormal->end();
+ const int numCellsNormal = cellsNormal->size();
+
+ // Count number of remaining cells (censored cells).
+ const int numCellsCensored = cells->size() - numCellsNormal;
+
+ // Get number of normal vertices.
+ assert(!mesh->getFactory().isNull());
+ Obj<mesh_type::numbering_type> vNumbering = mesh->getFactory()->getNumbering(mesh, "censored depth", 0);
+ assert(!vNumbering.isNull());
+ const int numVerticesNormal = vNumbering->size();
+
+ // Count number of remaining vertices (censored vertices).
+ const int numVerticesCensored = vertices->size() - numVerticesNormal;
+
+ _cellsNormal = ALE::Interval<point_type>(0, numCellsNormal);
+ _verticesNormal = ALE::Interval<point_type>(numCellsNormal, numCellsNormal+numVerticesNormal);
+ _verticesCensored = ALE::Interval<point_type>(numCellsNormal+numVerticesNormal, numCellsNormal+numVerticesNormal+numVerticesCensored);
+ _cellsCensored = ALE::Interval<point_type>(numCellsNormal+numVerticesNormal+numVerticesCensored,
+ numCellsNormal+numVerticesNormal+numVerticesCensored+numCellsCensored);
+ } else {
+ const int numCells = cells->size();
+ const int numVertices = vertices->size();
+ const int numEntities = numCells + numVertices;
+
+ _cellsNormal = ALE::Interval<point_type>(0, numCells);
+ _verticesNormal = ALE::Interval<point_type>(numCells, numCells+numVertices);
+ _cellsCensored = ALE::Interval<point_type>(numEntities, numEntities);
+ _verticesCensored = ALE::Interval<point_type>(numEntities, numEntities);
+ } // if/else
} // initialize
// ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.cc 2010-09-25 00:58:45 UTC (rev 17219)
@@ -87,7 +87,6 @@
ALE::ISieveVisitor::PointRetriever<mesh_type::sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
// Count number of cells.
- const int oldNumCells = cells->size();
int newNumCells = 0;
for (mesh_type::label_sequence::iterator c_iter = cells->begin(); c_iter != cellsEnd; ++c_iter) {
newNumCells += refiner.numNewCells(*c_iter);
@@ -184,6 +183,7 @@
_stratify(newMesh);
_calcNewOverlap(newMesh, mesh);
+ _createIntSections(newMesh, mesh, refiner);
} // _refine
// ----------------------------------------------------------------------
@@ -193,6 +193,214 @@
const Obj<mesh_type>& mesh,
CellRefinerTri3& refiner)
{ // _refineCensored
+ typedef Interval<point_type>::const_iterator interval_type;
+
+ assert(_orderOldMesh);
+ assert(_orderNewMesh);
+
+ // Calculate order in old mesh.
+ _orderOldMesh->initialize(mesh);
+
+ assert(!mesh.isNull());
+ assert(!newMesh.isNull());
+
+ // Get original mesh stuff.
+ const Obj<mesh_type::label_sequence>& cells = mesh->heightStratum(0);
+ assert(!cells.isNull());
+ const mesh_type::label_sequence::iterator cellsEnd = cells->end();
+
+ const Obj<mesh_type::label_sequence>& vertices = mesh->depthStratum(0);
+ assert(!vertices.isNull());
+
+ const Obj<mesh_type::sieve_type>& sieve = mesh->getSieve();
+ assert(!sieve.isNull());
+ ALE::ISieveVisitor::PointRetriever<mesh_type::sieve_type> cV(std::max(1, sieve->getMaxConeSize()));
+
+ // Count number of cells in censored depth (normal cells).
+ const Obj<mesh_type::label_sequence>& cellsNormal = mesh->getLabelStratum("censored depth", mesh->depth());
+ assert(!cellsNormal.isNull());
+ const mesh_type::label_sequence::iterator cellsNormalEnd = cellsNormal->end();
+ int newNumCellsNormal = 0;
+ for(mesh_type::label_sequence::iterator c_iter = cellsNormal->begin(); c_iter != cellsNormalEnd; ++c_iter)
+ newNumCellsNormal += refiner.numNewCells(*c_iter);
+
+ // Count number of remaining cells (other cells).
+ const int numSkip = cellsNormal->size();
+ mesh_type::label_sequence::iterator c_iter = cells->begin();
+ for (int i=0; i < numSkip; ++i)
+ ++c_iter;
+ int newNumCellsCensored = 0;
+ for (; c_iter != cellsEnd; ++c_iter)
+ newNumCellsCensored += refiner.numNewCells(*c_iter);
+
+ // Count number of normal vertices.
+ assert(!mesh->getFactory().isNull());
+ Obj<mesh_type::numbering_type> vNumbering = mesh->getFactory()->getNumbering(mesh, "censored depth", 0);
+ assert(!vNumbering.isNull());
+ const int oldNumVerticesNormal = vNumbering->size();
+
+ int counterBegin = newNumCellsNormal + oldNumVerticesNormal;
+ point_type curNewVertex = counterBegin;
+ for(mesh_type::label_sequence::iterator c_iter = cellsNormal->begin(); c_iter != cellsNormalEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ refiner.splitCell(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
+ } // for
+ const int newNumVerticesNormal = curNewVertex - counterBegin + oldNumVerticesNormal;
+
+ // Count number of remaining vertices (other vertices).
+ const int oldNumVerticesCensored = vertices->size() - oldNumVerticesNormal;
+ counterBegin = newNumCellsNormal + newNumVerticesNormal + oldNumVerticesCensored;
+ curNewVertex = counterBegin;
+ c_iter = cells->begin();
+ for (int i=0; i < numSkip; ++i)
+ ++c_iter;
+ for (; c_iter != cellsEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ refiner.splitCell(*c_iter, cV.getPoints(), cV.getSize(), &curNewVertex);
+ } // for
+ const int newNumVerticesCensored = curNewVertex - counterBegin + oldNumVerticesCensored;
+
+ _orderNewMesh->cellsNormal(0, newNumCellsNormal);
+ _orderNewMesh->verticesNormal(newNumCellsNormal, newNumCellsNormal+newNumVerticesNormal);
+ _orderNewMesh->verticesCensored(newNumCellsNormal+newNumVerticesNormal, newNumCellsNormal+newNumVerticesNormal+newNumVerticesCensored);
+ _orderNewMesh->cellsCensored(newNumCellsNormal+newNumVerticesNormal+newNumVerticesCensored,
+ newNumCellsNormal+newNumVerticesNormal+newNumVerticesCensored+newNumCellsCensored);
+
+ // Allocate chart for new sieve.
+ const Obj<mesh_type::sieve_type>& newSieve = newMesh->getSieve();
+ assert(!newSieve.isNull());
+ newSieve->setChart(mesh_type::sieve_type::chart_type(0, _orderNewMesh->cellsCensored().max()));
+
+ // Create new sieve with correct sizes for refined cells
+ point_type curNewCell = _orderNewMesh->cellsNormal().min();
+ interval_type::const_iterator oldCellsEnd = _orderOldMesh->cellsNormal().end();
+ for (interval_type::const_iterator c_iter=_orderOldMesh->cellsNormal().begin(); c_iter != oldCellsEnd; ++c_iter) {
+ // Set new cone and support sizes
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ const point_type* cone = cV.getPoints();
+ const int coneSize = cV.getSize();
+
+ const point_type* newCells;
+ int numNewCells = 0;
+ refiner.getNewCells(&newCells, &numNewCells, *c_iter, cone, coneSize, *_orderOldMesh, *_orderNewMesh);
+
+ for(int iCell=0; iCell < numNewCells; ++iCell, ++curNewCell) {
+ newSieve->setConeSize(curNewCell, coneSize);
+ for(int iVertex=0; iVertex < coneSize; ++iVertex) {
+ newSieve->addSupportSize(newCells[iCell*coneSize+iVertex], 1);
+ } // for
+ } // for
+ } // for
+ // Reset current new cell value and loop over censored cells.
+ curNewCell = _orderNewMesh->cellsCensored().min();
+ oldCellsEnd = _orderOldMesh->cellsCensored().end();
+ for (interval_type::const_iterator c_iter=_orderOldMesh->cellsCensored().begin(); c_iter != oldCellsEnd; ++c_iter) {
+ // Set new cone and support sizes
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ const point_type* cone = cV.getPoints();
+ const int coneSize = cV.getSize();
+
+ const point_type* newCells;
+ int numNewCells = 0;
+ refiner.getNewCells(&newCells, &numNewCells, *c_iter, cone, coneSize, *_orderOldMesh, *_orderNewMesh);
+
+ for(int iCell=0; iCell < numNewCells; ++iCell, ++curNewCell) {
+ newSieve->setConeSize(curNewCell, coneSize);
+ for(int iVertex=0; iVertex < coneSize; ++iVertex) {
+ newSieve->addSupportSize(newCells[iCell*coneSize+iVertex], 1);
+ } // for
+ } // for
+ } // for
+ newSieve->allocate();
+
+ // Create refined cells in new sieve.
+ curNewCell = _orderNewMesh->cellsNormal().min();
+ oldCellsEnd = _orderOldMesh->cellsNormal().end();
+ for (interval_type::const_iterator c_iter=_orderOldMesh->cellsNormal().begin(); c_iter != oldCellsEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ const point_type *cone = cV.getPoints();
+ const int coneSize = cV.getSize();
+
+ const point_type* newCells;
+ int numNewCells = 0;
+ refiner.getNewCells(&newCells, &numNewCells, *c_iter, cone, coneSize, *_orderOldMesh, *_orderNewMesh);
+
+ for(int iCell=0; iCell < numNewCells; ++iCell, ++curNewCell) {
+ newSieve->setCone(&newCells[iCell*coneSize], curNewCell);
+ } // for
+ } // for
+ curNewCell = _orderNewMesh->cellsCensored().min();
+ oldCellsEnd = _orderOldMesh->cellsCensored().end();
+ for (interval_type::const_iterator c_iter=_orderOldMesh->cellsCensored().begin(); c_iter != oldCellsEnd; ++c_iter) {
+ cV.clear();
+ sieve->cone(*c_iter, cV);
+ const point_type *cone = cV.getPoints();
+ const int coneSize = cV.getSize();
+
+ const point_type* newCells;
+ int numNewCells = 0;
+ refiner.getNewCells(&newCells, &numNewCells, *c_iter, cone, coneSize, *_orderOldMesh, *_orderNewMesh);
+
+ for(int iCell=0; iCell < numNewCells; ++iCell, ++curNewCell) {
+ newSieve->setCone(&newCells[iCell*coneSize], curNewCell);
+ } // for
+ } // for
+ newSieve->symmetrize();
+
+ // Set coordinates in refined mesh.
+ const Obj<mesh_type::real_section_type>& coordinates = mesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ const Obj<mesh_type::real_section_type>& newCoordinates = newMesh->getRealSection("coordinates");
+ assert(!newCoordinates.isNull());
+
+ const mesh_type::label_sequence::const_iterator verticesEnd = vertices->end();
+ assert(vertices->size() > 0);
+ const int spaceDim = coordinates->getFiberDimension(*vertices->begin());
+ assert(spaceDim > 0);
+ newCoordinates->setChart(mesh_type::sieve_type::chart_type(_orderNewMesh->verticesNormal().min(), _orderNewMesh->verticesCensored().max()));
+
+ const interval_type::const_iterator newVerticesEnd = _orderNewMesh->verticesCensored().end();
+ for (interval_type::const_iterator v_iter=_orderNewMesh->verticesNormal().begin(); v_iter != newVerticesEnd; ++v_iter) {
+ newCoordinates->setFiberDimension(*v_iter, spaceDim);
+ } // for
+ newCoordinates->allocatePoint();
+
+ interval_type::const_iterator oldVerticesEnd = _orderOldMesh->verticesNormal().end();
+ for (interval_type::const_iterator vOld_iter=_orderOldMesh->verticesNormal().begin(), vNew_iter=_orderNewMesh->verticesNormal().begin(); vOld_iter != oldVerticesEnd; ++vOld_iter, ++vNew_iter) {
+ std::cout << "Copy coordinates from old vertex " << *vOld_iter << " to new vertex " << *vNew_iter << std::endl;
+ newCoordinates->updatePoint(*vNew_iter, coordinates->restrictPoint(*vOld_iter));
+ } // for
+ oldVerticesEnd = _orderOldMesh->verticesCensored().end();
+ for (interval_type::const_iterator vOld_iter=_orderOldMesh->verticesCensored().begin(), vNew_iter=_orderNewMesh->verticesCensored().begin(); vOld_iter != oldVerticesEnd; ++vOld_iter, ++vNew_iter) {
+ std::cout << "Copy coordinates from old vertex " << *vOld_iter << " to new vertex " << *vNew_iter << std::endl;
+ newCoordinates->updatePoint(*vNew_iter, coordinates->restrictPoint(*vOld_iter));
+ } // for
+
+ refiner.setCoordsNewVertices(newCoordinates, coordinates);
+
+ _stratify(newMesh);
+ _calcNewOverlap(newMesh, mesh);
+ _createIntSections(newMesh, mesh, refiner);
+
+ // Create sensored depth
+ const ALE::Obj<Mesh::label_type>& censoredLabel = newMesh->createLabel("censored depth");
+ assert(!censoredLabel.isNull());
+
+ mesh_type::DepthVisitor depthVisitor(*newSieve, _orderNewMesh->verticesCensored().min(), *censoredLabel);
+
+ newSieve->roots(depthVisitor);
+ while(depthVisitor.isModified()) {
+ // FIX: Avoid the copy here somehow by fixing the traversal
+ std::vector<mesh_type::point_type> modifiedPoints(depthVisitor.getModifiedPoints().begin(), depthVisitor.getModifiedPoints().end());
+
+ depthVisitor.clear();
+ newSieve->support(modifiedPoints, depthVisitor);
+ } // while
} // _refineCensored
// ----------------------------------------------------------------------
@@ -238,6 +446,75 @@
} // _stratify
// ----------------------------------------------------------------------
+// Create integer sections in new mesh.
+void
+ALE::MeshRefiner::_createIntSections(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner)
+{ // _createIntSections
+ const ALE::Obj<std::set<std::string> >& sectionNames =
+ mesh->getIntSections();
+ const std::set<std::string>::const_iterator namesBegin =
+ sectionNames->begin();
+ const std::set<std::string>::const_iterator namesEnd =
+ sectionNames->end();
+ for (std::set<std::string>::const_iterator name=namesBegin;
+ name != namesEnd;
+ ++name) {
+ const ALE::Obj<mesh_type::int_section_type>& group = mesh->getIntSection(*name);
+ const ALE::Obj<mesh_type::int_section_type>& newGroup = newMesh->getIntSection(*name);
+ const mesh_type::int_section_type::chart_type& chart = group->getChart();
+
+ // :WARNING: Only implemented for integer sections containing vertices.
+
+ const point_type oldVerticesStart = _orderOldMesh->verticesNormal().min();
+ const point_type oldVerticesCensoredStart = _orderOldMesh->verticesCensored().min();
+ const point_type oldVerticesEnd = _orderOldMesh->verticesCensored().max();
+
+ const point_type newVerticesStart = _orderNewMesh->verticesNormal().min();
+ const point_type newVerticesCensoredStart = _orderNewMesh->verticesCensored().min();
+ const point_type newVerticesEnd = _orderNewMesh->verticesCensored().max();
+
+ newGroup->setChart(mesh_type::int_section_type::chart_type(newVerticesStart, newVerticesEnd));
+ const mesh_type::int_section_type::chart_type& newChart = newGroup->getChart();
+
+ std::cout << "VERTICES start: " << newVerticesStart << ", end: " << newVerticesEnd << std::endl;
+
+ const point_type chartMax = chart.max();
+ for (point_type pOld=chart.min(); pOld < chartMax; ++pOld) {
+ if (group->getFiberDimension(pOld)) {
+ if (_orderOldMesh->verticesNormal().hasPoint(pOld)) {
+ const point_type pNew = newVerticesStart + pOld - oldVerticesStart;
+ newGroup->setFiberDimension(pNew, 1);
+ } else if (_orderOldMesh->verticesCensored().hasPoint(pOld)) {
+ const point_type pNew = newVerticesCensoredStart + pOld - oldVerticesCensoredStart;
+ newGroup->setFiberDimension(pNew, 1);
+ } else {
+ throw ALE::Exception("Creating integer sections during refinement containing entities other than vertices not implemented.");
+ } // if/else
+ } // if
+ } // for
+ refiner.groupAddNewVertices(newGroup, group);
+
+ newGroup->allocatePoint();
+ for (point_type pOld=chart.min(); pOld < chartMax; ++pOld) {
+ if (group->getFiberDimension(pOld)) {
+ if (_orderOldMesh->verticesNormal().hasPoint(pOld)) {
+ const point_type pNew = newVerticesStart + pOld - oldVerticesStart;
+ newGroup->updatePoint(pNew, group->restrictPoint(pOld));
+ } else if (_orderOldMesh->verticesCensored().hasPoint(pOld)) {
+ const point_type pNew = newVerticesCensoredStart + pOld - oldVerticesCensoredStart;
+ newGroup->updatePoint(pNew, group->restrictPoint(pOld));
+ } else {
+ throw ALE::Exception("Creating integer sections during refinement containing entities other than vertices not implemented.");
+ } // if/else
+ } // if
+ } // for
+ refiner.groupSetNewVertices(newGroup, group);
+ } // for
+} // _createIntSections
+
+// ----------------------------------------------------------------------
// Calculate new overlap.
void
ALE::MeshRefiner::_calcNewOverlap(const Obj<mesh_type>& newMesh,
Modified: short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/libsrc/topology/MeshRefiner.hh 2010-09-25 00:58:45 UTC (rev 17219)
@@ -94,6 +94,17 @@
void _calcNewOverlap(const Obj<mesh_type>& newMesh,
const Obj<mesh_type>& mesh);
+ /** Create integer sections in new mesh.
+ *
+ * :WARNING: Only implemented for integer sections containing vertices.
+ *
+ * @param newMesh New (refined) mesh.
+ * @param mesh Current (unrefined) mesh with integer sections.
+ */
+ void _createIntSections(const Obj<mesh_type>& newMesh,
+ const Obj<mesh_type>& mesh,
+ CellRefinerTri3& refiner);
+
// PRIVATE MEMBERS //////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineUniform.cc 2010-09-25 00:58:45 UTC (rev 17219)
@@ -88,7 +88,6 @@
const ALE::Obj<SieveMesh::label_type>& newMaterials =
newSieveMesh->createLabel("material-id");
-
for (SieveMesh::label_sequence::const_iterator c_iter = cellsBegin,
cNew_iter = newCellsBegin;
c_iter != cellsEnd;
@@ -100,66 +99,6 @@
newSieveMesh->setValue(newMaterials, *cNew_iter, material);
} // for
- // Recreate groups, assuming vertex groups
- const int numNewVertices = newSieveMesh->depthStratum(0)->size();
- const int numNewCells = newSieveMesh->heightStratum(0)->size();
- const ALE::Obj<std::set<std::string> >& sectionNames =
- sieveMesh->getIntSections();
-
-#if 0
- ALE::MeshBuilder<SieveMesh>::CellRefiner<SieveMesh,edge_type>::edge_map_type& edge2vertex =
- refiner.getEdgeToVertex();
-
- const std::set<std::string>::const_iterator namesBegin =
- sectionNames->begin();
- const std::set<std::string>::const_iterator namesEnd =
- sectionNames->end();
- for (std::set<std::string>::const_iterator name=namesBegin;
- name != namesEnd;
- ++name) {
- const ALE::Obj<Mesh::IntSection>& group = sieveMesh->getIntSection(*name);
- const ALE::Obj<Mesh::IntSection>& newGroup =
- newSieveMesh->getIntSection(*name);
- const Mesh::IntSection::chart_type& chart = group->getChart();
-
- newGroup->setChart(Mesh::IntSection::chart_type(numNewCells,
- numNewCells + numNewVertices));
- const Mesh::IntSection::chart_type& newChart = newGroup->getChart();
-
- const int chartMax = chart.max();
- for (int p = chart.min(), pNew = newChart.min(); p < chartMax; ++p, ++pNew) {
- if (group->getFiberDimension(p))
- newGroup->setFiberDimension(pNew, 1);
- } // for
- const std::map<edge_type, point_type>::const_iterator edge2VertexEnd =
- edge2vertex.end();
- for (std::map<edge_type, point_type>::const_iterator e_iter=edge2vertex.begin();
- e_iter != edge2VertexEnd;
- ++e_iter) {
- const point_type vertexA = e_iter->first.first;
- const point_type vertexB = e_iter->first.second;
- if (group->getFiberDimension(vertexA) && group->getFiberDimension(vertexB))
- if (group->restrictPoint(vertexA)[0] == group->restrictPoint(vertexB)[0])
- newGroup->setFiberDimension(e_iter->second, 1);
- } // for
-
- newGroup->allocatePoint();
- for (int p=chart.min(), pNew = newChart.min(); p < chartMax; ++p, ++pNew) {
- if (group->getFiberDimension(p))
- newGroup->updatePoint(pNew, group->restrictPoint(p));
- } // for
- for (std::map<edge_type, point_type>::const_iterator e_iter=edge2vertex.begin();
- e_iter != edge2VertexEnd;
- ++e_iter) {
- const point_type vertexA = e_iter->first.first;
- const point_type vertexB = e_iter->first.second;
-
- if (group->getFiberDimension(vertexA) && group->getFiberDimension(vertexB))
- if (group->restrictPoint(vertexA)[0] == group->restrictPoint(vertexB)[0])
- newGroup->updatePoint(e_iter->second, group->restrictPoint(vertexA));
- } // for
- } // for
-#endif
} // refine
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2010-09-25 00:58:45 UTC (rev 17219)
@@ -261,7 +261,7 @@
const int numPoints = groupField->size();
int_array points(numPoints);
int i = 0;
- const int offset = ("vertex" == groupType) ? numCells : 0;
+ const int offset = ("vertex" == groupType) ? data.numCells : 0;
for(chart_type::const_iterator c_iter = chart.begin();
c_iter != chart.end();
++c_iter)
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.hh 2010-09-25 00:58:45 UTC (rev 17219)
@@ -52,8 +52,8 @@
CPPUNIT_TEST( testConstructor );
CPPUNIT_TEST( testRefineTri3Level2 );
+ CPPUNIT_TEST( testRefineTri3Level2Fault1 );
#if 0
- CPPUNIT_TEST( testRefineTri3Level2Fault1 );
CPPUNIT_TEST( testRefineTet4Level2 );
CPPUNIT_TEST( testRefineTet4Level2Fault1 );
#endif
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/data/MeshDataCohesiveTri3Level2Fault1.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/data/MeshDataCohesiveTri3Level2Fault1.cc 2010-09-24 22:46:00 UTC (rev 17218)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/data/MeshDataCohesiveTri3Level2Fault1.cc 2010-09-25 00:58:45 UTC (rev 17219)
@@ -49,9 +49,6 @@
0.0, -1.0,
0.0, 0.0,
0.0, 1.0,
- 0.0, -1.0,
- 0.0, 0.0,
- 0.0, 1.0,
-0.5, -0.5,
0.0, -0.5,
-0.5, 0.0,
@@ -62,35 +59,36 @@
0.5, 0.0,
0.5, 0.5,
0.0, 0.5,
+ 0.0, -1.0,
+ 0.0, 0.0,
+ 0.0, 1.0,
0.0, -0.5,
0.0, 0.5,
};
const int pylith::topology::MeshDataCohesiveTri3Level2Fault1::_cells[] = {
- 0, 11, 13,
- 11, 12, 13,
- 1, 12, 11,
- 2, 13, 12,
- 2, 14, 13,
- 14, 15, 13,
- 3, 15, 14,
- 0, 13, 15,
- 6, 16, 18,
- 16, 17, 18,
- 5, 17, 16,
- 4, 18, 17,
- 6, 18, 20,
- 18, 19, 20,
- 4, 19, 18,
- 7, 20, 19,
-
+ 0, 8, 10,
+ 8, 9, 10,
+ 1, 9, 8,
+ 2, 10, 9,
+ 2, 11, 10,
+ 11, 12, 10,
+ 3, 12, 11,
+ 0, 10, 12,
+ 6, 13, 15,
+ 13, 14, 15,
+ 5, 14, 13,
+ 4, 15, 14,
+ 6, 15, 17,
+ 15, 16, 17,
+ 4, 16, 15,
+ 7, 17, 16,
};
const int pylith::topology::MeshDataCohesiveTri3Level2Fault1::_cellsCohesive[] = {
- 1, 12, 5, 16, 8, 21,
- 12, 2, 16, 6, 21, 9,
- 2, 14, 6, 20, 9, 22,
- 14, 3, 20, 7, 22, 10,
-
+ 1, 9, 5, 13, 18, 21,
+ 9, 2, 13, 6, 21, 19,
+ 2, 11, 6, 17, 19, 22,
+ 11, 3, 17, 7, 22, 20,
};
const int pylith::topology::MeshDataCohesiveTri3Level2Fault1::_materialIds[] = {
1, 1, 1, 1, 1, 1, 1, 1,
@@ -105,10 +103,10 @@
};
const int pylith::topology::MeshDataCohesiveTri3Level2Fault1::_groups[] = {
- 0, 1, 3, 5, 7, 11, 15,
- 1, 4, 5, 17,
+ 0, 1, 3, 5, 7, 8, 12,
+ 1, 4, 5, 14,
0, 4,
- 1, 2, 3, 5, 6, 7, 8, 9, 10, 12, 14, 16, 20, 21, 22,
+ 1, 2, 3, 5, 6, 7, 9, 11, 13, 17, 18, 19, 20, 21, 22,
};
const char* pylith::topology::MeshDataCohesiveTri3Level2Fault1::_groupNames[] = {
More information about the CIG-COMMITS
mailing list