00001 #include "SundanceRivaraElement.hpp"
00002 #include "SundanceRivaraEdge.hpp"
00003 #include "SundanceRivaraNode.hpp"
00004 #include "SundanceRivaraMesh.hpp"
00005 #include "SundanceOut.hpp"
00006
00007
00008
00009 using namespace Sundance::Rivara;
00010 using namespace Teuchos;
00011 using std::endl;
00012
00013
00014 Element::Element(RivaraMesh* mesh,
00015 const RCP<Node>& a,
00016 const RCP<Node>& b,
00017 const RCP<Node>& c,
00018 int ownerProc,
00019 int label)
00020 : label_(label),
00021 nodes_(tuple(a,b,c)),
00022 edges_(3),
00023 faces_(0),
00024 edgeSigns_(3),
00025 ownerProc_(ownerProc)
00026 {
00027 for (int i=0; i<3; i++)
00028 {
00029 edges_[i] = mesh->tryEdge(nodes_[i], nodes_[(i+1)%3], edgeSigns_[i]);
00030 edges_[i]->addConnectingElement(this);
00031 nodes_[i]->addConnectingElement(this);
00032
00033 nodes_[i]->addConnectingEdge(edges_[i].get());
00034 nodes_[(i+1) % 3]->addConnectingEdge(edges_[i].get());
00035 }
00036 }
00037
00038 Element::Element(RivaraMesh* mesh,
00039 const RCP<Node>& a,
00040 const RCP<Node>& b,
00041 const RCP<Node>& c,
00042 const RCP<Node>& d,
00043 int ownerProc,
00044 int label)
00045 : label_(label),
00046 nodes_(tuple(a,b,c,d)),
00047 edges_(6),
00048 faces_(4),
00049 edgeSigns_(6),
00050 ownerProc_(ownerProc)
00051 {
00052
00053 for (int i=0; i<4; i++)
00054 {
00055 nodes_[i]->addConnectingElement(this);
00056 }
00057
00058 int k=0;
00059 for (int i=0; i<4; i++)
00060 {
00061 for (int j=0; j<i; j++)
00062 {
00063 edges_[k] = mesh->tryEdge(nodes_[i], nodes_[j], edgeSigns_[k]);
00064 edges_[k]->addConnectingElement(this);
00065 nodes_[i]->addConnectingEdge(edges_[k].get());
00066 nodes_[j]->addConnectingEdge(edges_[k].get());
00067 k++;
00068 }
00069 }
00070
00071 faces_[0] = mesh->tryFace(a,b,d);
00072 faces_[1] = mesh->tryFace(b,c,d);
00073 faces_[2] = mesh->tryFace(a,d,c);
00074 faces_[3] = mesh->tryFace(c,b,a);
00075 }
00076
00077 int Element::longestEdgeIndex() const
00078 {
00079 int e = 0;
00080 double L = -1.0;
00081 for (int i=0; i<edges_.length(); i++)
00082 {
00083 if (edges_[i]->length() > L) {e = i; L = edges_[i]->length();}
00084 }
00085 return e;
00086 }
00087
00088 bool Element::hasHangingNode() const
00089 {
00090 bool hasHangingNode = false;
00091 for (int i=0; i<edges_.length(); i++)
00092 {
00093 if (edges_[i]->hasChildren()) hasHangingNode = true;
00094 }
00095 return hasHangingNode;
00096 }
00097
00098 void Element::refine(RivaraMesh* mesh, double maxArea)
00099 {
00100
00101 if (hasChildren())
00102 {
00103 dynamic_cast<Element*>(left())->refine(mesh, maxArea);
00104 dynamic_cast<Element*>(right())->refine(mesh, maxArea);
00105 return;
00106 }
00107
00108
00109
00110
00111 if (!hasHangingNode() && volume() < maxArea)
00112 {
00113 return;
00114 }
00115 if (!hasHangingNode() && maxArea < 0.0)
00116 {
00117 return;
00118 }
00119
00120
00121
00122
00123
00124 int e = longestEdgeIndex();
00125
00126 Element* sub1;
00127 Element* sub2;
00128
00129 if (nodes_.length()==3)
00130 {
00131
00132
00133
00134
00135 RCP<Node> a, b, c;
00136 if (edgeSigns_[e] > 0)
00137 {
00138 a = edges_[e]->node(0);
00139 b = edges_[e]->node(1);
00140 }
00141 else
00142 {
00143 b = edges_[e]->node(0);
00144 a = edges_[e]->node(1);
00145 }
00146
00147 for (int i=0; i<nodes_.length(); i++)
00148 {
00149 if (nodes_[i].get() == a.get() || nodes_[i].get() == b.get()) continue;
00150 c = nodes_[i];
00151 break;
00152 }
00153
00154
00155
00156 RCP<Node> mid = edges_[e]->bisect(mesh);
00157
00158
00159 sub1 = new Element(mesh, a, mid, c, ownerProc_, label_);
00160 sub2 = new Element(mesh, c, mid, b, ownerProc_, label_);
00161 }
00162 else
00163 {
00164
00165
00166
00167
00168
00169 RCP<Node> a, b, c, d;
00170 if (edgeSigns_[e] > 0)
00171 {
00172 a = edges_[e]->node(0);
00173 b = edges_[e]->node(1);
00174 }
00175 else
00176 {
00177 b = edges_[e]->node(0);
00178 a = edges_[e]->node(1);
00179 }
00180
00181 for (int i=0; i<edges_.length(); i++)
00182 {
00183 const RCP<Node>& node1 = edges_[i]->node(0);
00184 const RCP<Node>& node2 = edges_[i]->node(1);
00185 if (node1.get()==a.get() || node1.get()==b.get()
00186 || node2.get()==a.get() || node2.get()==b.get())
00187 {
00188 continue;
00189 }
00190 if (edgeSigns_[i] > 0)
00191 {
00192 c = edges_[i]->node(0);
00193 d = edges_[i]->node(1);
00194 }
00195 else
00196 {
00197 d = edges_[i]->node(0);
00198 c = edges_[i]->node(1);
00199 }
00200 break;
00201 }
00202
00203
00204
00205 RCP<Node> mid = edges_[e]->bisect(mesh);
00206
00207
00208 sub1 = new Element(mesh, a, mid, c, d, ownerProc_, label_);
00209 sub2 = new Element(mesh, b, c, mid, d, ownerProc_, label_);
00210
00211
00212
00213
00214 if (label_ != -1)
00215 {
00216 RCP<Face> abc = mesh->getFace(a,b,c);
00217 RCP<Face> abd = mesh->getFace(a,b,d);
00218 RCP<Face> acd = mesh->getFace(a,c,d);
00219 RCP<Face> bcd = mesh->getFace(b,c,d);
00220
00221 RCP<Face> amd = mesh->getFace(a,mid,d);
00222 RCP<Face> amc = mesh->getFace(a,mid,c);
00223 RCP<Face> bmc = mesh->getFace(b,mid,c);
00224 RCP<Face> bmd = mesh->getFace(b,mid,d);
00225
00226 amd->setLabel(abd->label());
00227 amc->setLabel(abc->label());
00228
00229 bmd->setLabel(abd->label());
00230 bmc->setLabel(abc->label());
00231 }
00232
00233 }
00234
00235 sub1->setParent(this);
00236 sub2->setParent(this);
00237
00238 setChildren(sub1, sub2);
00239
00240
00241
00242 Array<Element*> others;
00243 edges_[e]->getUnrefinedCofacets(others);
00244 for (int i=0; i<others.length(); i++)
00245 {
00246 Element* other = others[i];
00247 if (other==0) continue;
00248 mesh->refinementSet().push(other);
00249 mesh->refinementAreas().push(-1.0);
00250 }
00251
00252 if (sub1->hasHangingNode() || maxArea > 0.0)
00253 sub1->refine(mesh, maxArea);
00254 if (sub2->hasHangingNode() || maxArea > 0.0)
00255 sub2->refine(mesh, maxArea);
00256 }
00257
00258 namespace Sundance
00259 {
00260 inline Point cross3(const Point& a, const Point& b)
00261 {
00262 return Point(a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2],
00263 a[0]*b[1]-a[1]*b[0]);
00264 }
00265
00266 inline double cross2(const Point& a, const Point& b)
00267 {
00268 return a[0]*b[1] - a[1]*b[0];
00269 }
00270 }
00271
00272 double Element::volume() const
00273 {
00274 if (nodes_.length()==3)
00275 {
00276 return 0.5*::fabs(cross2(nodes_[2]->pt()-nodes_[0]->pt(),
00277 nodes_[1]->pt()-nodes_[0]->pt()));
00278 }
00279 else
00280 {
00281 Point AB = nodes_[1]->pt() - nodes_[0]->pt();
00282 Point AC = nodes_[2]->pt() - nodes_[0]->pt();
00283 Point AD = nodes_[3]->pt() - nodes_[0]->pt();
00284 return 1.0/6.0 * fabs( AB * cross3(AC, AD) );
00285 }
00286 }
00287
00288
00289 Array<int> Element::showNodes() const
00290 {
00291 Array<int> rtn(nodes_.length());
00292 for (int i=0; i<nodes_.length(); i++) rtn[i]=nodes_[i]->localIndex();
00293 return rtn;
00294 }
00295
00296 bool Element::hasNoEdgeLabels() const
00297 {
00298 for (int i=0; i<3; i++)
00299 {
00300 if (edges_[i]->label() != -1) return false;
00301 }
00302 return true;
00303 }
00304
00305 bool Element::hasNoFaceLabels() const
00306 {
00307 for (int i=0; i<4; i++)
00308 {
00309 if (faces_[i]->label() != -1) return false;
00310 }
00311 return true;
00312 }