00001 #include "SundanceExtrusionMeshTransformation.hpp"
00002 #include "PlayaExceptions.hpp"
00003 #include "SundanceOut.hpp"
00004
00005 using namespace Sundance;
00006 using namespace Sundance;
00007
00008 using namespace Teuchos;
00009 using namespace Sundance;
00010
00011
00012
00013 Mesh ExtrusionMeshTransformation::apply(const Mesh& inputMesh) const
00014 {
00015 TEUCHOS_TEST_FOR_EXCEPTION(inputMesh.spatialDim() != 2, std::runtime_error,
00016 "ExtrusionMeshTransformation::applyLocal() given mesh with "
00017 "dimension " << inputMesh.spatialDim() << ". The "
00018 "extrusion filter expects a 2D mesh as input");
00019
00020
00021 Mesh outputMesh = createMesh(3, inputMesh.comm());
00022
00023
00024
00025
00026 int nPts = inputMesh.numCells(0);
00027
00028
00029 int pointCount = 0;
00030 outputMesh.estimateNumVertices(nzLevels_ * nPts);
00031
00032 for (int i=0; i<nPts; i++)
00033 {
00034 const Point& p = inputMesh.nodePosition(i);
00035 for (int j=0; j<=nzLevels_; j++, pointCount++)
00036 {
00037 double z = z0_ + (z1_-z0_)*j/((double) nzLevels_);
00038 outputMesh.addVertex(pointCount, Point(p[0], p[1], z), 0, 0);
00039 }
00040 }
00041
00042
00043
00044 int nTri = inputMesh.numCells(2);
00045
00046
00047
00048
00049 outputMesh.estimateNumElements(3*nTri*nzLevels_);
00050
00051 Array<int> facetNodes(3);
00052 Array<int> facetOrientations(3);
00053
00054 int tetCount = 0;
00055
00056 TEUCHOS_TEST_FOR_EXCEPTION(inputMesh.cellType(2) != TriangleCell,
00057 std::runtime_error,
00058 "ExtrusionMeshTransformation::applyLocal() detected a "
00059 "non-triangular mesh");
00060 for (int i=0; i<nTri; i++)
00061 {
00062 inputMesh.getFacetArray(2, i, 0, facetNodes, facetOrientations);
00063
00064 SUNDANCE_OUT(this->verb()==3,
00065 "triangle=" << i << " facet nodes are " << facetNodes);
00066
00067
00068
00069 std::sort(facetNodes.begin(), facetNodes.end());
00070
00071 SUNDANCE_OUT(this->verb()==3,
00072 "triangle=" << i << " sorted facet nodes are "
00073 << facetNodes);
00074
00075 int a0 = facetNodes[0]*(nzLevels_+1);
00076 int b0 = facetNodes[1]*(nzLevels_+1);
00077 int c0 = facetNodes[2]*(nzLevels_+1);
00078
00079 for (int j=0; j<nzLevels_; j++, tetCount+=3)
00080 {
00081 int a = a0 + j;
00082 int b = b0 + j;
00083 int c = c0 + j;
00084 int a1 = a+1;
00085 int b1 = b+1;
00086 int c1 = c+1;
00087 outputMesh.addElement(tetCount, tuple(a, a1, b1, c1), 0, 0);
00088 outputMesh.addElement(tetCount+1, tuple(a, b, b1, c1), 0, 0);
00089 outputMesh.addElement(tetCount+2, tuple(a, b, c, c1), 0, 0);
00090 }
00091 }
00092
00093 return outputMesh;
00094
00095 }
00096