00001 #include "SundanceRivaraMesh.hpp" 00002 #include "SundanceRivaraEdge.hpp" 00003 #include "SundanceRivaraElement.hpp" 00004 #include "SundanceRivaraNode.hpp" 00005 #include "Teuchos_StrUtils.hpp" 00006 #include <fstream> 00007 #include "Teuchos_TimeMonitor.hpp" 00008 #include "SundanceOut.hpp" 00009 00010 using namespace Sundance::Rivara; 00011 using namespace Sundance; 00012 using Sundance::Map; 00013 00014 using namespace Teuchos; 00015 using std::endl; 00016 00017 00018 00019 00020 static Time& refKernelTimer() 00021 { 00022 static RCP<Time> rtn 00023 = TimeMonitor::getNewTimer("mesh refinement kernel"); 00024 return *rtn; 00025 } 00026 00027 RivaraMesh::RivaraMesh(int dim, const MPIComm& comm) 00028 : spatialDim_(dim), nextGID_(0), 00029 nodes_(), edges_(), elements_(), nodeToEdgeMap_() 00030 {;} 00031 00032 void RivaraMesh::refine() 00033 { 00034 TimeMonitor timer(refKernelTimer()); 00035 while (refinementSet().size() > 0) 00036 { 00037 Element* next = refinementSet().top(); 00038 refinementSet().pop(); 00039 double maxArea = refinementAreas().top(); 00040 refinementAreas().pop(); 00041 next->refine(this, maxArea); 00042 } 00043 } 00044 00045 int RivaraMesh::addNode(const RCP<Node>& node) 00046 { 00047 int lid = nodes_.length(); 00048 node->setLocalIndex(lid); 00049 nodes_.append(node); 00050 nodeToEdgeMap_.append(Map<int, int>()); 00051 nextGID()++; 00052 return lid; 00053 } 00054 00055 00056 int RivaraMesh::addVertex( 00057 int globalIndex, const Point& x, 00058 int ownerProcID, int label) 00059 { 00060 RCP<Node> node = rcp(new Node(globalIndex, x, ownerProcID, label)); 00061 return addNode(node); 00062 } 00063 00064 void RivaraMesh::addElement(const RCP<Element>& tri) 00065 { 00066 elements_.append(tri); 00067 } 00068 00069 00070 int RivaraMesh::addElement( 00071 int globalIndex, 00072 const Array<int>& vertexGIDs, 00073 int ownerProc, 00074 int label) 00075 { 00076 int lid = elements_.size(); 00077 RCP<Element> elem; 00078 00079 switch(vertexGIDs.size()) 00080 { 00081 case 3: 00082 elem = rcp(new Element(this, 00083 nodes_[vertexGIDs[0]], 00084 nodes_[vertexGIDs[1]], 00085 nodes_[vertexGIDs[2]], 00086 ownerProc,label)); 00087 break; 00088 case 4: 00089 elem = rcp(new Element(this, 00090 nodes_[vertexGIDs[0]], 00091 nodes_[vertexGIDs[1]], 00092 nodes_[vertexGIDs[2]], 00093 nodes_[vertexGIDs[3]], 00094 ownerProc,label)); 00095 break; 00096 default: 00097 TEUCHOS_TEST_FOR_EXCEPT(1); 00098 } 00099 elements_.append(elem); 00100 return lid; 00101 00102 } 00103 00104 RCP<Edge> RivaraMesh::tryEdge(const RCP<Node>& a, 00105 const RCP<Node>& b, 00106 int& edgeSign) 00107 { 00108 int i = a->localIndex(); 00109 int j = b->localIndex(); 00110 00111 if (nodeToEdgeMap_[i].containsKey(j)) 00112 { 00113 int k = nodeToEdgeMap_[i].get(j); 00114 edgeSign = 1; 00115 return edges_[k]; 00116 } 00117 else if (nodeToEdgeMap_[j].containsKey(i)) 00118 { 00119 int k = nodeToEdgeMap_[j].get(i); 00120 edgeSign = -1; 00121 return edges_[k]; 00122 } 00123 else 00124 { 00125 RCP<Edge> rtn = rcp(new Edge(a,b)); 00126 edgeSign = 1; 00127 int k = edges_.length(); 00128 edges_.append(rtn); 00129 nodeToEdgeMap_[i].put(j, k); 00130 return rtn; 00131 } 00132 } 00133 00134 00135 RCP<Face> RivaraMesh::tryFace( 00136 const RCP<Node>& a, 00137 const RCP<Node>& b, 00138 const RCP<Node>& c) 00139 { 00140 FaceNodes f(a,b,c); 00141 00142 int faceLID; 00143 if (faceToLIDMap_.containsKey(f)) 00144 { 00145 faceLID = faceToLIDMap_[f]; 00146 } 00147 else 00148 { 00149 faceLID = faces_.size(); 00150 RCP<Face> newFace = rcp(new Face(a,b,c)); 00151 faceToLIDMap_.put(newFace->nodes(), faceLID); 00152 faces_.append(newFace); 00153 } 00154 00155 return faces_[faceLID]; 00156 } 00157 00158 00159 00160 const RCP<Face>& RivaraMesh::getFace( 00161 const RCP<Node>& a, 00162 const RCP<Node>& b, 00163 const RCP<Node>& c) const 00164 { 00165 FaceNodes f(a,b,c); 00166 return faces_[faceToLIDMap_.get(f)]; 00167 } 00168 00169 00170 00171 int RivaraMesh::numElements() const 00172 { 00173 int numLeaves = 0; 00174 for (int i=0; i<elements_.length(); i++) 00175 { 00176 numLeaves += elements_[i]->numLeaves(); 00177 } 00178 return numLeaves; 00179 } 00180 00181 int RivaraMesh::spatialDim() const 00182 { 00183 return spatialDim_; 00184 } 00185 00186