[cig-commits] r17268 - short/3D/PyLith/trunk/libsrc/topology
brad at geodynamics.org
brad at geodynamics.org
Wed Oct 13 22:06:24 PDT 2010
Author: brad
Date: 2010-10-13 22:06:23 -0700 (Wed, 13 Oct 2010)
New Revision: 17268
Added:
short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.cc
short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.hh
short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.cc
short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.hh
Log:
Started work on uniform refinement for quad4 cells.
Added: short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.cc 2010-10-14 05:06:23 UTC (rev 17268)
@@ -0,0 +1,423 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "CellRefinerQuad4.hh" // implementation of class methods
+
+#include "MeshOrder.hh" // USES MeshOrder
+
+#include <cassert> // USES assert()
+
+#include <iostream> // TEMPORARY
+// ----------------------------------------------------------------------
+// Constructor
+ALE::CellRefinerQuad4::CellRefinerQuad4(const mesh_type& mesh) :
+ RefineFace4Edges2(mesh)
+{ // constructor
+ assert(2 == mesh.getDimension());
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+ALE::CellRefinerQuad4::~CellRefinerQuad4(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get number of refined cells for each original cell.
+int
+ALE::CellRefinerQuad4::numNewCells(const point_type cell)
+{ // numNewCells
+ switch (_cellType(cell)) {
+ case QUADRILATERAL:
+ return 4;
+ case LINE_COHESIVE_LAGRANGE:
+ return 2;
+ default:
+ assert(0);
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+} // numNewCells
+
+// ----------------------------------------------------------------------
+// Split cell into smaller cells of same type.
+void
+ALE::CellRefinerQuad4::splitCell(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex)
+{ // splitCell
+ assert(curNewVertex);
+
+ int numEdges = 0;
+ const EdgeType* edges;
+
+ int numFaces = 0;
+ const FaceType* faces;
+
+ switch (_cellType(cell)) {
+ case QUADRILATERAL:
+ _edges_QUADRILATERAL(&edges, &numEdges, cone, coneSize);
+ _faces_QUADRILATERAL(&faces, &numFaces, cone, coneSize);
+ break;
+ case LINE_COHESIVE_LAGRANGE:
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize);
+ break;
+ default:
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+
+ for(int iEdge=0; iEdge < numEdges; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ // if vertex does not exist
+ std::cout << "Edge: " << edges[iEdge] << ", new vertex: " << *curNewVertex << std::endl;
+ _edgeToVertex[edges[iEdge]] = *curNewVertex;
+ ++(*curNewVertex);
+ } // if
+ } // for
+
+ for(int iFace=0; iFace < numFaces; ++iFace) {
+ if (_faceToVertex.find(faces[iFace]) == _faceToVertex.end()) {
+ // if vertex does not exist
+ std::cout << "Face: " << faces[iFace] << ", new vertex: " << *curNewVertex << std::endl;
+ _faceToVertex[faces[iFace]] = *curNewVertex;
+ ++(*curNewVertex);
+ } // if
+ } // for
+} // splitCell
+
+// ----------------------------------------------------------------------
+// Split cell into smaller cells of same type.
+void
+ALE::CellRefinerQuad4::splitCellUncensored(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex)
+{ // splitCellUncensored
+ assert(curNewVertex);
+
+ int numEdges = 0;
+ const EdgeType* edges;
+
+ const bool uncensored = true;
+
+ switch (_cellType(cell)) {
+ case QUADRILATERAL:
+ // No censored vertices on normal cell.
+ break;
+ case LINE_COHESIVE_LAGRANGE:
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize, uncensored);
+ break;
+ default:
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+
+ for(int iEdge=0; iEdge < numEdges; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ // if vertex does not exist
+ std::cout << "Edge: " << edges[iEdge] << ", new vertex: " << *curNewVertex << std::endl;
+ _edgeToVertex[edges[iEdge]] = *curNewVertex;
+ ++(*curNewVertex);
+ } // if
+ } // for
+
+} // splitCellUncensored
+
+// ----------------------------------------------------------------------
+// Get refined cells.
+void
+ALE::CellRefinerQuad4::getNewCells(const point_type** cells,
+ int* numCells,
+ const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ const MeshOrder& orderOldMesh,
+ const MeshOrder& orderNewMesh)
+{ // getNewCells
+ assert(cells);
+ assert(numCells);
+
+ switch (_cellType(cell)) {
+ case QUADRILATERAL: {
+ const int coneVertexOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
+ _newCells_QUADRILATERAL(cells, numCells, cone, coneSize, coneVertexOffset);
+ break;
+ } // QUADRILATERAL
+ case LINE_COHESIVE_LAGRANGE: {
+ 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:
+ throw ALE::Exception("Unknown cell type.");
+ } // switch
+} // getNewCells
+
+// ----------------------------------------------------------------------
+// Get cell type.
+ALE::CellRefinerQuad4::CellEnum
+ALE::CellRefinerQuad4::_cellType(const point_type cell)
+{ // _cellType
+ assert(!_mesh.getSieve().isNull());
+
+ switch (_mesh.getSieve()->getConeSize(cell)) {
+ case 4:
+ return QUADRILATERAL;
+ case 6:
+ return LINE_COHESIVE_LAGRANGE;
+ case 0: {
+ std::ostringstream msg;
+ std::cerr << "Internal error. Cone size for mesh point " << cell << " is zero. May be a vertex.";
+ assert(0);
+ throw ALE::Exception("Could not determine cell type during uniform global refinement.");
+ } // case 0
+ default : {
+ std::ostringstream msg;
+ std::cerr << "Internal error. Unknown cone size for mesh point " << cell << ". Unknown cell type.";
+ assert(0);
+ throw ALE::Exception("Could not determine cell type during uniform global refinement.");
+ } // default
+ } // switch
+} // _cellType
+
+// ----------------------------------------------------------------------
+// Get edges of quadrilateral cell.
+void
+ALE::CellRefinerQuad4::_edges_QUADRILATERAL(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize)
+{ // _edges_QUADRILATERAL
+ static EdgeType quadEdges[4];
+
+ assert(coneSize == 4);
+ quadEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ quadEdges[1] = EdgeType(std::min(cone[1], cone[2]), std::max(cone[1], cone[2]));
+ quadEdges[2] = EdgeType(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
+ quadEdges[3] = EdgeType(std::min(cone[3], cone[0]), std::max(cone[3], cone[0]));
+ *numEdges = 4;
+ *edges = quadEdges;
+} // _edges_QUADRILATERAL
+
+// ----------------------------------------------------------------------
+// Get edges of line cohesive cell with Lagrange multipler vertices.
+void
+ALE::CellRefinerQuad4::_edges_LINE_COHESIVE_LAGRANGE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize,
+ const bool uncensored)
+{ // _edges_LINE_COHESIVE_LAGRANGE
+ if (uncensored) {
+ // Include all edges
+ static EdgeType lineEdges[3];
+
+ assert(coneSize == 6);
+ lineEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ lineEdges[1] = EdgeType(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
+ lineEdges[2] = EdgeType(std::min(cone[4], cone[5]), std::max(cone[4], cone[5]));
+ *numEdges = 3;
+ *edges = lineEdges;
+ } else {
+ // Omit edges with censored (Lagrange multiplier) vertices.
+ static EdgeType lineEdges[2];
+
+ assert(coneSize == 6);
+ lineEdges[0] = EdgeType(std::min(cone[0], cone[1]), std::max(cone[0], cone[1]));
+ lineEdges[1] = EdgeType(std::min(cone[2], cone[3]), std::max(cone[2], cone[3]));
+ *numEdges = 2;
+ *edges = lineEdges;
+ } // if/else
+} // _edges_LINE_COHESIVE_LAGRANGE
+
+// ----------------------------------------------------------------------
+// Get faces of quadrilateral cell.
+void
+ALE::CellRefinerQuad4::_faces_QUADRILATERAL(const FaceType** faces,
+ int* numFaces,
+ const point_type cone[],
+ const int coneSize)
+{ // _faces_QUADRILATERAL
+ static FaceType quadFaces[1];
+
+ assert(coneSize == 4);
+
+ int sortedCone[4];
+ for (int i=0; i < 4; ++i)
+ sortedCone[i] = cone[i];
+ std::sort(sortedCone, sortedCone+coneSize);
+ const point_type pMin = sortedCone[0];
+
+ if (pMin == cone[0]) {
+ if (cone[1] < cone[3]) {
+ quadFaces[0] = FaceType(cone[0], cone[1], cone[2], cone[3]);
+ } else {
+ quadFaces[0] = FaceType(cone[0], cone[3], cone[2], cone[1]);
+ } // if/else
+
+ } else if (pMin == cone[1]) {
+ if (cone[2] < cone[0]) {
+ quadFaces[0] = FaceType(cone[1], cone[2], cone[3], cone[0]);
+ } else {
+ quadFaces[0] = FaceType(cone[1], cone[0], cone[3], cone[2]);
+ } // if/else
+
+ } else if (pMin == cone[2]) {
+ if (cone[3] < cone[1]) {
+ quadFaces[0] = FaceType(cone[2], cone[3], cone[0], cone[1]);
+ } else {
+ quadFaces[0] = FaceType(cone[2], cone[1], cone[3], cone[0]);
+ } // if/else
+
+ } else if (pMin == cone[3]) {
+ if (cone[0] < cone[2]) {
+ quadFaces[0] = FaceType(cone[3], cone[0], cone[1], cone[2]);
+ } else {
+ quadFaces[0] = FaceType(cone[3], cone[2], cone[1], cone[0]);
+ } // if/else
+ } else {
+ assert(0);
+ throw ALE::Exception("Could not determine quad face orientation during uniform global refinement.");
+ } // if/else
+ *numFaces = 1;
+ *faces = quadFaces;
+} // _faces_QUADRILATERAL
+
+// ----------------------------------------------------------------------
+// Get new cells from refinement of a triangular cell.
+void
+ ALE::CellRefinerQuad4::_newCells_QUADRILATERAL(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset)
+{ // _newCells_QUADRILATERAL
+ const int coneSizeQuad4 = 4;
+ const int numEdgesQuad4 = 4;
+ const int numFacesQuad4 = 1;
+ const int numNewCells = 4;
+ const int numNewVertices = 5;
+
+ int numEdges = 0;
+ const EdgeType *edges;
+ _edges_QUADRILATERAL(&edges, &numEdges, cone, coneSize);
+ assert(numEdgesQuad4 == numEdges);
+
+ int numFaces = 0;
+ const FaceType *faces;
+ _faces_QUADRILATERAL(&faces, &numFaces, cone, coneSize);
+ assert(numFacesQuad4 == numFaces);
+
+ static point_type quadCells[numNewCells*coneSizeQuad4];
+ point_type newVertices[numNewVertices];
+ int iNewVertex = 0;
+ for(int iEdge=0; iEdge < numEdgesQuad4; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ throw ALE::Exception("Missing edge in refined mesh");
+ } // if
+ newVertices[iNewVertex++] = _edgeToVertex[edges[iEdge]];
+ } // for
+ for(int iFace=0; iFace < numFacesQuad4; ++iFace) {
+ if (_faceToVertex.find(faces[iFace]) == _faceToVertex.end()) {
+ throw ALE::Exception("Missing face in refined mesh");
+ } // if
+ newVertices[iNewVertex++] = _faceToVertex[faces[iFace]];
+ } // for
+
+ // new cell 0
+ quadCells[0*4+0] = cone[0] + coneVertexOffset;
+ quadCells[0*4+1] = newVertices[0];
+ quadCells[0*4+2] = newVertices[4];
+ quadCells[0*4+3] = newVertices[3];
+
+ // new cell 1
+ quadCells[1*4+0] = cone[1] + coneVertexOffset;
+ quadCells[1*4+1] = newVertices[1];
+ quadCells[1*4+2] = newVertices[4];
+ quadCells[1*4+3] = newVertices[0];
+
+ // new cell 2
+ quadCells[2*4+0] = cone[3] + coneVertexOffset;
+ quadCells[2*4+1] = newVertices[3];
+ quadCells[2*4+2] = newVertices[4];
+ quadCells[2*4+3] = newVertices[2];
+
+ // new cell 3
+ quadCells[3*4+0] = cone[2] + coneVertexOffset;
+ quadCells[3*4+1] = newVertices[3];
+ quadCells[3*4+2] = newVertices[4];
+ quadCells[3*4+3] = newVertices[1];
+
+ *numCells = numNewCells;
+ *cells = quadCells;
+} // _newCells_QUADRILATERAL
+
+// ----------------------------------------------------------------------
+// Get new cells from refinement of a line cohseive cell with Lagrange
+// multiplier vertices.
+void
+ALE::CellRefinerQuad4::_newCells_LINE_COHESIVE_LAGRANGE(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffsetNormal,
+ const int coneVertexOffsetCensored)
+{ // _newCells_LINE_COHESIVE_LAGRANGE
+ const int coneSizeLine6 = 6;
+ const int numEdgesLine6 = 3;
+ const int numNewCells = 2;
+ const int numNewVertices = 3;
+
+ int numEdges = 0;
+ const EdgeType *edges;
+ _edges_LINE_COHESIVE_LAGRANGE(&edges, &numEdges, cone, coneSize, true);
+ assert(numEdgesLine6 == numEdges);
+
+ static point_type lineCells[numNewCells*coneSizeLine6];
+ point_type newVertices[numNewVertices];
+ for(int iEdge=0, iNewVertex=0; iEdge < numEdgesLine6; ++iEdge) {
+ if (_edgeToVertex.find(edges[iEdge]) == _edgeToVertex.end()) {
+ throw ALE::Exception("Missing edge in refined mesh");
+ } // if
+ newVertices[iNewVertex++] = _edgeToVertex[edges[iEdge]];
+ } // for
+
+ // new cell 0
+ lineCells[0*6+0] = cone[0] + coneVertexOffsetNormal;
+ lineCells[0*6+1] = newVertices[0];
+ lineCells[0*6+2] = cone[2] + coneVertexOffsetNormal;
+ lineCells[0*6+3] = newVertices[1];
+ 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] + coneVertexOffsetNormal;
+ lineCells[1*6+2] = newVertices[1];
+ lineCells[1*6+3] = cone[3] + coneVertexOffsetNormal;
+ lineCells[1*6+4] = newVertices[2];
+ lineCells[1*6+5] = cone[5] + coneVertexOffsetCensored;
+
+ *numCells = 2;
+ *cells = lineCells;
+} // _newCells_LINE_COHESIVE_LAGRANGE
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/CellRefinerQuad4.hh 2010-10-14 05:06:23 UTC (rev 17268)
@@ -0,0 +1,210 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/CellRefinerQuad4.hh
+ *
+ * @brief Object for quad4 refinement of cells.
+ */
+
+#if !defined(pylith_topology_cellrefinerquad4_hh)
+#define pylith_topology_cellrefinerquad4_hh
+
+// Include directives ---------------------------------------------------
+#include "RefineFace4Edges2.hh" // ISA RefineFace4Edges2
+
+// CellRefinerQuad4 ------------------------------------------------------
+/// Object for quad4 refinement of cells.
+class ALE::CellRefinerQuad4 : public RefineFace4Edges2
+{ // CellRefinerQuad4
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Constructor
+ *
+ * @param mesh Finite-element mesh.
+ */
+ CellRefinerQuad4(const mesh_type& mesh);
+
+ /// Destructor
+ ~CellRefinerQuad4(void);
+
+ /** Get number of refined cells for each original cell.
+ *
+ * @param cell Original cell.
+ *
+ * @returns Number of refined cells.
+ */
+ int numNewCells(const point_type cell);
+
+ /** Split cell into smaller cells of same type. Do not create
+ * censored vertices on censored cells.
+ *
+ * @param cell Original cell.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param curNewVertex Value for next new vertex.
+ */
+ void splitCell(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex);
+
+ /** Split cell into smaller cells of same type. Create only censored
+ * vertices on censored cells.
+ *
+ * @param cell Original cell.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param curNewVertex Value for next new vertex.
+ */
+ void splitCellUncensored(const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ point_type* curNewVertex);
+
+ /** Get refined cells.
+ *
+ * @param cells Vertices in refined cells (refined mesh).
+ * @param numCells Number of refined cells.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param orderOldMesh Order in old mesh.
+ * @param orderNewMesh Order in new mesh.
+ */
+ void getNewCells(const point_type** cells,
+ int* numCells,
+ const point_type cell,
+ const point_type cone[],
+ const int coneSize,
+ const MeshOrder& orderOldMesh,
+ const MeshOrder& orderNewMesh);
+
+// PRIVATE ENUMS ////////////////////////////////////////////////////////
+private :
+
+ enum CellEnum {
+ QUADRILATERAL, // Normal quadrilateral cell
+ LINE_COHESIVE_LAGRANGE, // Cohesive cell with Lagrange multiplier vertices
+ };
+
+// PRIVATE METHODS //////////////////////////////////////////////////////
+private :
+
+ /** Get cell type.
+ *
+ * @param cell Cell in original mesh.
+ * @returns Cell type.
+ */
+ CellEnum _cellType(const point_type cell);
+
+ /** Get edges of quadrilateral cell.
+ *
+ * @param edges Edges of cell.
+ * @param numEdges Number of edges.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ */
+ void _edges_QUADRILATERAL(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize);
+
+ /** Get edges of line cohesive cell with Lagrange multipler vertices.
+ *
+ * @param edges Edges of cell.
+ * @param numEdges Number of edges.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param uncensored True if including edges with censored vertices.
+ */
+ void _edges_LINE_COHESIVE_LAGRANGE(const EdgeType** edges,
+ int* numEdges,
+ const point_type cone[],
+ const int coneSize,
+ const bool uncensored =false);
+
+ /** Get faces of quadrilateral cell.
+ *
+ * @param faces Faces of cell.
+ * @param numFaces Number of faces.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ */
+ void _faces_QUADRILATERAL(const FaceType** faces,
+ int* numFaces,
+ const point_type cone[],
+ const int coneSize);
+
+ /** Get faces of line cohesive cell with Lagrange multipler vertices.
+ *
+ * @param faces Faces of cell.
+ * @param numFaces Number of faces.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @param uncensored True if including faces with censored vertices.
+ */
+ void _faces_LINE_COHESIVE_LAGRANGE(const FaceType** faces,
+ int* numFaces,
+ const point_type cone[],
+ const int coneSize,
+ const bool uncensored =false);
+
+ /** Get new cells from refinement of a quadrilateral cell.
+ *
+ * @param cells Vertices in refined cells (refined mesh).
+ * @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.
+ */
+ void _newCells_QUADRILATERAL(const point_type** cells,
+ int *numCells,
+ const point_type cone[],
+ const int coneSize,
+ const int coneVertexOffset);
+
+ /** Get new cells from refinement of a line cohseive cell with
+ * Lagrange multiplier vertices.
+ *
+ * @param cells Vertices in refined cells (refined mesh).
+ * @param numCells Number of refined cells.
+ * @param cone Vertices in cell (original mesh).
+ * @param coneSize Number of vertices in cell.
+ * @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 coneVertexOffsetNormal,
+ const int coneVertexOffsetCensored);
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ CellRefinerQuad4(const CellRefinerQuad4&); ///< Not implemented
+ const CellRefinerQuad4& operator=(const CellRefinerQuad4&); ///< Not implemented
+
+}; // CellRefinerQuad4
+
+#endif // pylith_topology_cellrefinerquad4_hh
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.cc (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.cc 2010-10-14 05:06:23 UTC (rev 17268)
@@ -0,0 +1,276 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "RefineFace4Edges2.hh" // implementation of class methods
+
+#include "MeshOrder.hh" // USES MeshOrder
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Constructor
+ALE::RefineFace4Edges2::RefineFace4Edges2(const mesh_type& mesh) :
+ _mesh(mesh)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+ALE::RefineFace4Edges2::~RefineFace4Edges2(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set coordinates of new vertices.
+void
+ALE::RefineFace4Edges2::setCoordsNewVertices(const ALE::Obj<mesh_type::real_section_type>& newCoordsSection,
+ const ALE::Obj<mesh_type::real_section_type>& oldCoordsSection)
+{ // setCoordsNewVertices
+ assert(!newCoordsSection.isNull());
+ assert(!oldCoordsSection.isNull());
+
+ double coordinatesVertex[3];
+
+ assert(_edgeToVertex.size() > 0);
+ const int spaceDim = newCoordsSection->getFiberDimension(_edgeToVertex.begin()->second);
+ assert(spaceDim > 0 && spaceDim <= 3);
+
+ 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;
+
+ assert(spaceDim == oldCoordsSection->getFiberDimension(edgeVertexA));
+ assert(spaceDim == oldCoordsSection->getFiberDimension(edgeVertexB));
+ assert(spaceDim == newCoordsSection->getFiberDimension(newVertex));
+
+ const mesh_type::real_section_type::value_type* coordsA = oldCoordsSection->restrictPoint(edgeVertexA);
+ const mesh_type::real_section_type::value_type* coordsB = oldCoordsSection->restrictPoint(edgeVertexB);
+ for (int i=0; i < spaceDim; ++i)
+ coordinatesVertex[i] = 0.5*(coordsA[i] + coordsB[i]);
+
+ newCoordsSection->updatePoint(newVertex, coordinatesVertex);
+ } // for
+} // setCoordsNewVertices
+
+// ----------------------------------------------------------------------
+// Add space for new vertices in group.
+void
+ALE::RefineFace4Edges2::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::RefineFace4Edges2::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
+
+// ----------------------------------------------------------------------
+// Add new vertices to label.
+void
+ALE::RefineFace4Edges2::labelAddNewVertices(const ALE::Obj<mesh_type>& newMesh,
+ const ALE::Obj<mesh_type>& oldMesh,
+ const char* labelName)
+{ // labelAddNewVertices
+ assert(!newMesh.isNull());
+ assert(!oldMesh.isNull());
+
+ const Obj<mesh_type::label_sequence>& oldLabelVertices = oldMesh->getLabelStratum(labelName, 0);
+ assert(!oldLabelVertices.isNull());
+
+ const Obj<mesh_type::label_type>& oldLabel = oldMesh->getLabel(labelName);
+ assert(!oldLabel.isNull());
+ const Obj<mesh_type::label_type>& newLabel = newMesh->getLabel(labelName);
+ assert(!newLabel.isNull());
+
+ const int defaultValue = -999;
+
+ 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;
+
+ const int valueA = oldMesh->getValue(oldLabel, edgeVertexA, defaultValue);
+ const int valueB = oldMesh->getValue(oldLabel, edgeVertexB, defaultValue);
+
+ if (valueA != defaultValue && valueA == valueB) {
+ newMesh->setValue(newLabel, newVertex, valueA);
+ } // if
+ } // for
+} // labelAddNewVertices
+
+// ----------------------------------------------------------------------
+// Calculate new overlap.
+void
+ALE::RefineFace4Edges2::overlapAddNewVertices(const Obj<mesh_type>& newMesh,
+ const MeshOrder& orderNewMesh,
+ const Obj<mesh_type>& oldMesh,
+ const MeshOrder& orderOldMesh)
+{ // overlapAddNewVertices
+ assert(!newMesh.isNull());
+ assert(!oldMesh.isNull());
+
+ Obj<mesh_type::send_overlap_type> newSendOverlap = newMesh->getSendOverlap();
+ assert(!newSendOverlap.isNull());
+ Obj<mesh_type::recv_overlap_type> newRecvOverlap = newMesh->getRecvOverlap();
+ assert(!newRecvOverlap.isNull());
+ const Obj<mesh_type::send_overlap_type>& oldSendOverlap = oldMesh->getSendOverlap();
+ assert(!oldSendOverlap.isNull());
+ const Obj<mesh_type::recv_overlap_type>& oldRecvOverlap = oldMesh->getRecvOverlap();
+ assert(!oldRecvOverlap.isNull());
+
+
+ // Check edges in edgeToVertex for both endpoints sent to same process
+ // Put it in section with point being the lowest numbered vertex and value (other endpoint, new vertex)
+ Obj<ALE::Section<point_type, EdgeType> > newVerticesSection = new ALE::Section<point_type, EdgeType>(oldMesh->comm());
+ assert(!newVerticesSection.isNull());
+ std::map<EdgeType, std::vector<int> > bndryEdgeToRank;
+
+ const int localOffset = orderNewMesh.verticesNormal().min() - orderOldMesh.verticesNormal().min();
+
+ for(std::map<EdgeType, point_type>::const_iterator e_iter = _edgeToVertex.begin(); e_iter != _edgeToVertex.end(); ++e_iter) {
+ const point_type left = e_iter->first.first;
+ const point_type right = e_iter->first.second;
+
+ if (oldSendOverlap->capContains(left) && oldSendOverlap->capContains(right)) {
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanksSeq = oldSendOverlap->support(left);
+ assert(!leftRanksSeq.isNull());
+ std::list<int> leftRanks(leftRanksSeq->begin(), leftRanksSeq->end());
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = oldSendOverlap->support(right);
+ assert(!rightRanks.isNull());
+ std::list<int> ranks;
+ std::set_intersection(leftRanks.begin(), leftRanks.end(), rightRanks->begin(), rightRanks->end(),
+ std::insert_iterator<std::list<int> >(ranks, ranks.begin()));
+
+ if(ranks.size()) {
+ newVerticesSection->addFiberDimension(std::min(e_iter->first.first, e_iter->first.second)+localOffset, 1);
+ for(std::list<int>::const_iterator r_iter = ranks.begin(); r_iter != ranks.end(); ++r_iter) {
+ bndryEdgeToRank[e_iter->first].push_back(*r_iter);
+ } // for
+ } // if
+ } // if
+ } // for
+ newVerticesSection->allocatePoint();
+ const ALE::Section<point_type, EdgeType>::chart_type& chart = newVerticesSection->getChart();
+
+ for(ALE::Section<point_type, EdgeType>::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+ typedef ALE::Section<point_type, EdgeType>::value_type value_type;
+ const point_type p = *c_iter;
+ const int dim = newVerticesSection->getFiberDimension(p);
+ int v = 0;
+ value_type* values = (dim > 0) ? new value_type[dim] : 0;
+
+ for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end() && v < dim; ++e_iter) {
+ if (std::min(e_iter->first.first, e_iter->first.second)+localOffset == p) {
+ values[v++] = EdgeType(std::max(e_iter->first.first, e_iter->first.second)+localOffset, _edgeToVertex[e_iter->first]);
+ } // if
+ } // for
+ newVerticesSection->updatePoint(p, values);
+ delete [] values;
+ } // for
+ // Copy across overlap
+ typedef ALE::Pair<int, point_type> overlap_point_type;
+ Obj<ALE::Section<overlap_point_type, EdgeType> > overlapVertices = new ALE::Section<overlap_point_type, EdgeType>(oldMesh->comm());
+
+ ALE::Pullback::SimpleCopy::copy(newSendOverlap, newRecvOverlap, newVerticesSection, overlapVertices);
+ // Merge by translating edge to local points, finding edge in _edgeToVertex, and adding (local new vetex, remote new vertex) to overlap
+ for(std::map<EdgeType, std::vector<int> >::const_iterator e_iter = bndryEdgeToRank.begin(); e_iter != bndryEdgeToRank.end(); ++e_iter) {
+ const point_type localPoint = _edgeToVertex[e_iter->first];
+
+ for(std::vector<int>::const_iterator r_iter = e_iter->second.begin(); r_iter != e_iter->second.end(); ++r_iter) {
+ point_type remoteLeft = -1, remoteRight = -1;
+ const int rank = *r_iter;
+
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& leftRanks = newSendOverlap->support(e_iter->first.first+localOffset);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator lr_iter = leftRanks->begin(); lr_iter != leftRanks->end(); ++lr_iter) {
+ if (rank == *lr_iter) {
+ remoteLeft = lr_iter.color();
+ break;
+ } // if
+ } // for
+ const Obj<mesh_type::send_overlap_type::traits::supportSequence>& rightRanks = newSendOverlap->support(e_iter->first.second+localOffset);
+ for(mesh_type::send_overlap_type::traits::supportSequence::iterator rr_iter = rightRanks->begin(); rr_iter != rightRanks->end(); ++rr_iter) {
+ if (rank == *rr_iter) {
+ remoteRight = rr_iter.color();
+ break;
+ } // if
+ } // for
+ const point_type remoteMin = std::min(remoteLeft, remoteRight);
+ const point_type remoteMax = std::max(remoteLeft, remoteRight);
+ const int remoteSize = overlapVertices->getFiberDimension(overlap_point_type(rank, remoteMin));
+ const EdgeType *remoteVals = overlapVertices->restrictPoint(overlap_point_type(rank, remoteMin));
+ point_type remotePoint = -1;
+
+ for(int d = 0; d < remoteSize; ++d) {
+ if (remoteVals[d].first == remoteMax) {
+ remotePoint = remoteVals[d].second;
+ break;
+ } // if
+ } // for
+ newSendOverlap->addArrow(localPoint, rank, remotePoint);
+ newRecvOverlap->addArrow(rank, localPoint, remotePoint);
+ } // for
+ } // for
+
+ oldSendOverlap->view("OLD SEND OVERLAP");
+ oldRecvOverlap->view("OLD RECV OVERLAP");
+ newSendOverlap->view("NEW SEND OVERLAP");
+ newRecvOverlap->view("NEW RECV OVERLAP");
+} // overlapAddNewVertces
+
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/topology/RefineFace4Edges2.hh 2010-10-14 05:06:23 UTC (rev 17268)
@@ -0,0 +1,198 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard, U.S. Geological Survey
+// Charles A. Williams, GNS Science
+// Matthew G. Knepley, University of Chicago
+//
+// This code was developed as part of the Computational Infrastructure
+// for Geodynamics (http://geodynamics.org).
+//
+// Copyright (c) 2010 University of California, Davis
+//
+// See COPYING for license information.
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/RefineFace4Edges2.hh
+ *
+ * @brief Object for refinement of cells via refinement of edges
+ * comprised of two vertices.
+ */
+
+#if !defined(pylith_topology_refineface4edges2_hh)
+#define pylith_topology_refineface4edges2_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+#include <list> // USES std::pair
+
+// RefineFace4Edges2 ------------------------------------------------------
+/// Object for tri3 refinement of cells.
+class ALE::RefineFace4Edges2
+{ // RefineFace4Edges2
+protected:
+
+ typedef IMesh<> mesh_type;
+ typedef mesh_type::point_type point_type;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Constructor
+ *
+ * @param mesh Finite-element mesh.
+ */
+ RefineFace4Edges2(const mesh_type& mesh);
+
+ /// Destructor
+ ~RefineFace4Edges2(void);
+
+ /** Set coordinates of new vertices.
+ *
+ * @param newCoordsSection Coordinates of vertices in new mesh.
+ * @param oldCoordsSection Coordinates of vertices in original mesh.
+ */
+ 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);
+
+ /** Add new vertices to label.
+ *
+ * @param newMesh Mesh with refined cells.
+ * @param oldMesh Original mesh.
+ * @param labelName Name of label.
+ */
+ void labelAddNewVertices(const ALE::Obj<mesh_type>& newMesh,
+ const ALE::Obj<mesh_type>& oldMesh,
+ const char* labelName);
+
+ /** Calculate new overlap.
+ *
+ * @param newMesh New (refined) mesh.
+ * @param orderNewMesh Order in new mesh.
+ * @param oldMesh Current (unrefined) mesh with overlap.
+ * @param orderOldMesh Order in old mesh.
+ */
+ void overlapAddNewVertices(const Obj<mesh_type>& newMesh,
+ const MeshOrder& orderNewMesh,
+ const Obj<mesh_type>& oldMesh,
+ const MeshOrder& orderOldMesh);
+
+// PROTECTED TYPEDEFS ///////////////////////////////////////////////////
+protected :
+
+ template<typename Point>
+ class Edge : public std::pair<Point, Point> {
+ public:
+ Edge(void) : std::pair<Point, Point>() {};
+ Edge(const Point l) : std::pair<Point, Point>(l, l) {};
+ Edge(const Point l, const Point r) : std::pair<Point, Point>(l, r) {};
+ ~Edge(void) {};
+ friend std::ostream& operator<<(std::ostream& stream, const Edge& edge) {
+ stream << "(" << edge.first << ", " << edge.second << ")";
+ return stream;
+ };
+ };
+ typedef Edge<point_type> EdgeType;
+ typedef std::map<EdgeType, point_type> edge_map_type;
+
+ template<typename Point>
+ class Face {
+ public:
+ Face(void);
+ Face(const Point p) {
+ _points[0] = p;
+ _points[1] = p;
+ _points[2] = p;
+ _points[3] = p;
+ };
+ Face(const Point p0,
+ const Point p1,
+ const Point p2,
+ const Point p3) {
+ _points[0] = p0;
+ _points[1] = p1;
+ _points[2] = p2;
+ _points[3] = p3;
+ };
+ ~Face(void) {};
+ friend bool operator==(const Face& a, const Face& b) {
+ const bool result =
+ a._points[0] == b._points[0] &&
+ a._points[1] == b._points[1] &&
+ a._points[2] == b._points[2] &&
+ a._points[3] == b._points[3];
+ return result;
+ };
+ friend bool operator<(const Face& a, const Face& b) {
+ if (a._points[0] < b._points[0]) {
+ return true;
+ } else if (a._points[0] == b._points[0]) {
+ if (a._points[1] < b._points[1]) {
+ return true;
+ } else if (a._points[1] == b._points[1]) {
+ if (a._points[2] < b._points[2]) {
+ return true;
+ } else if (a._points[2] == b._points[2]) {
+ if (a._points[3] < b._points[3]) {
+ return true;
+ } // if
+ } // if/else
+ } // if/else
+ } // if/else
+
+ return false;
+ };
+ friend std::ostream& operator<<(std::ostream& stream, const Face& face) {
+ stream << "(" << face._points[0]
+ << ", " << face._points[1]
+ << ", " << face._points[2]
+ << ", " << face._points[3]
+ << ")";
+ return stream;
+ };
+ private:
+ Point _points[4];
+ };
+ typedef Face<point_type> FaceType;
+ typedef std::map<FaceType, point_type> face_map_type;
+
+// PROTECTED MEMBERS ////////////////////////////////////////////////////
+protected :
+
+ const mesh_type& _mesh;
+ edge_map_type _edgeToVertex;
+ face_map_type _faceToVertex;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ RefineFace4Edges2(const RefineFace4Edges2&); ///< Not implemented
+ const RefineFace4Edges2& operator=(const RefineFace4Edges2&); ///< Not implemented
+
+}; // RefineFace4Edges2
+
+#endif // pylith_topology_refineface4edges2_hh
+
+
+// End of file
More information about the CIG-COMMITS
mailing list