1 //
2 // triangle.cc
3 //
4 // Copyright (C) 1996 Limit Point Systems, Inc.
5 //
6 // Author: Curtis Janssen <cljanss@limitpt.com>
7 // Maintainer: LPS
8 //
9 // This file is part of the SC Toolkit.
10 //
11 // The SC Toolkit is free software; you can redistribute it and/or modify
12 // it under the terms of the GNU Library General Public License as published by
13 // the Free Software Foundation; either version 2, or (at your option)
14 // any later version.
15 //
16 // The SC Toolkit is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU Library General Public License for more details.
20 //
21 // You should have received a copy of the GNU Library General Public License
22 // along with the SC Toolkit; see the file COPYING.LIB. If not, write to
23 // the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
24 //
25 // The U.S. Government is granted a limited license as per AL 91-7.
26 //
27
28 #ifdef __GNUC__
29 #pragma implementation
30 #endif
31
32 #include <util/misc/formio.h>
33 #include <util/keyval/keyval.h>
34 #include <math/isosurf/triangle.h>
35 #include <math/scmat/vector3.h>
36
37 using namespace std;
38 using namespace sc;
39
40 /////////////////////////////////////////////////////////////////////////
41 // Triangle
42 ///////////////////////////////////////////////////////////////////////////
43 // Here is the layout of the vertices and edges in the triangles: |
44 // |
45 // vertex(1) (r=0, s=1) |
46 // + |
47 // / \ _edges[1].vertex(_orientation1) |
48 // / \ |
49 // / \ |
50 // / \ |
51 // / \ |
52 // / \ |
53 // _edges[0] / \ _edges[1] |
54 // (r = 0) / \ (1-r-s = 0) |
55 // / \ |
56 // / \ |
57 // / \ |
58 // / \ _edges[1].vertex(!_orientation1)|
59 // / \ |
60 // vertex(0)+---------------------------+ |
61 // (r=0, s=0) _edges[2] (s = 0) vertex(2) (r=1, s=0) |
62 // |
63 // Zienkiewicz and Taylor, "The Finite Element Method", 4th Ed, Vol 1, |
64 // use |
65 // L1 = 1 - r - s |
66 // L2 = r, |
67 // L3 = s. |
68 // I also use these below. |
69 ///////////////////////////////////////////////////////////////////////////
70
Triangle(const Ref<Edge> & v1,const Ref<Edge> & v2,const Ref<Edge> & v3,unsigned int orientation0)71 Triangle::Triangle(const Ref<Edge>& v1, const Ref<Edge>& v2, const Ref<Edge>& v3,
72 unsigned int orientation0)
73 {
74 _orientation0 = orientation0;
75
76 _edges[0] = v1;
77 _edges[1] = v2;
78 _edges[2] = v3;
79
80 // edge 0 corresponds to r = 0
81 // edge 1 corresponds to (1-r-s) = 0
82 // edge 2 corresponds to s = 0
83 // edge 0 vertex _orientation0 is (r=0,s=0)
84 // edge 1 vertex _orientation1 is (r=0,s=1)
85 // edge 2 vertex _orientation2 is (r=1,s=0)
86 // edge 0 vertex (1-_orientation0) is edge 1 vertex _orientation1
87 // edge 1 vertex (1-_orientation1) is edge 2 vertex _orientation2
88 // edge 2 vertex (1-_orientation2) is edge 0 vertex _orientation0
89
90 Ref<Edge> *e = _edges;
91
92 // swap edges 1 and 2 if necessary
93 if (e[0]->vertex(1-_orientation0) != e[1]->vertex(0)) {
94 if (e[0]->vertex(1-_orientation0) != e[1]->vertex(1)) {
95 e[1] = v3;
96 e[2] = v2;
97 }
98 }
99
100 // compute the orientation of _edge[1]
101 if (e[0]->vertex(1-_orientation0) == e[1]->vertex(0)) {
102 _orientation1 = 0;
103 }
104 else {
105 _orientation1 = 1;
106 }
107
108 // compute the orientation of _edge[2]
109 if (e[1]->vertex(1-_orientation1) == e[2]->vertex(0)) {
110 _orientation2 = 0;
111 }
112 else {
113 _orientation2 = 1;
114 }
115
116 // make sure that the triangle is really a triangle
117 if ( e[0]->vertex(1-_orientation0) != e[1]->vertex(_orientation1)
118 || e[1]->vertex(1-_orientation1) != e[2]->vertex(_orientation2)
119 || e[2]->vertex(1-_orientation2) != e[0]->vertex(_orientation0))
120 {
121 ExEnv::errn() << "Triangle: given edges that don't form a triangle" << endl;
122 abort();
123 }
124
125 _order = 1;
126 _vertices = new Ref<Vertex>[3];
127
128 _vertices[TriInterpCoef::ijk_to_index(_order, 0, 0)] = vertex(0);
129 _vertices[TriInterpCoef::ijk_to_index(0, 0, _order)] = vertex(1);
130 _vertices[TriInterpCoef::ijk_to_index(0, _order, 0)] = vertex(2);
131 }
132
~Triangle()133 Triangle::~Triangle()
134 {
135 if (_vertices) delete[] _vertices;
136 }
137
138 Ref<Vertex>
vertex(int i)139 Triangle::vertex(int i)
140 {
141 return _edges[i]->vertex(orientation(i));
142 }
143
144 unsigned int
orientation(const Ref<Edge> & e) const145 Triangle::orientation(const Ref<Edge>& e) const
146 {
147 if (e == _edges[0]) return orientation(0);
148 if (e == _edges[1]) return orientation(1);
149 return orientation(2);
150 }
151
152 int
contains(const Ref<Edge> & e) const153 Triangle::contains(const Ref<Edge>& e) const
154 {
155 if (_edges[0] == e) return 1;
156 if (_edges[1] == e) return 1;
157 if (_edges[2] == e) return 1;
158 return 0;
159 }
160
flat_area()161 double Triangle::flat_area()
162 {
163 double a = _edges[0]->straight_length();
164 double b = _edges[1]->straight_length();
165 double c = _edges[2]->straight_length();
166 double a2 = a*a;
167 double b2 = b*b;
168 double c2 = c*c;
169 return 0.25 * sqrt( 2.0 * (c2*b2 + c2*a2 + a2*b2)
170 - c2*c2 - b2*b2 - a2*a2);
171 }
172
add_vertices(std::set<Ref<Vertex>> & set)173 void Triangle::add_vertices(std::set<Ref<Vertex> >&set)
174 {
175 for (int i=0; i<3; i++) set.insert(_edges[i]->vertex(orientation(i)));
176 }
177
add_edges(std::set<Ref<Edge>> & set)178 void Triangle::add_edges(std::set<Ref<Edge> >&set)
179 {
180 for (int i=0; i<3; i++) set.insert(_edges[i]);
181 }
182
183 void
interpolate(double r,double s,const Ref<Vertex> & result,SCVector3 & dA)184 Triangle::interpolate(double r,double s,const Ref<Vertex>&result,SCVector3&dA)
185 {
186 TriInterpCoefKey key(_order, r, s);
187 Ref<TriInterpCoef> coef = new TriInterpCoef(key);
188 interpolate(coef, r, s, result, dA);
189 }
190
191 void
interpolate(const Ref<TriInterpCoef> & coef,double r,double s,const Ref<Vertex> & result,SCVector3 & dA)192 Triangle::interpolate(const Ref<TriInterpCoef>& coef,
193 double r, double s,
194 const Ref<Vertex>&result, SCVector3&dA)
195 {
196 unsigned int i, j, k;
197
198 //double L1 = 1 - r - s;
199 //double L2 = r;
200 //double L3 = s;
201
202 SCVector3 tmp(0.0);
203 SCVector3 x_s(0.0);
204 SCVector3 x_r(0.0);
205 for (i=0; i<=_order; i++) {
206 for (j=0; j <= _order-i; j++) {
207 k = _order - i - j;
208 int index = TriInterpCoef::ijk_to_index(i,j,k);
209 tmp += coef->coef(i,j,k)*_vertices[index]->point();
210 x_s += coef->sderiv(i,j,k)*_vertices[index]->point();
211 x_r += coef->rderiv(i,j,k)*_vertices[index]->point();
212 }
213 }
214 result->set_point(tmp);
215
216 if (result->has_normal()) {
217 for (i=0; i<_order; i++) {
218 for (j=0; j <= _order-i; j++) {
219 k = _order - i - j;
220 int index = TriInterpCoef::ijk_to_index(i,j,k);
221 tmp += coef->coef(i,j,k)*_vertices[index]->point();
222 }
223 }
224 result->set_normal(tmp);
225 }
226
227 // Find the surface element
228 dA = x_r.cross(x_s);
229 }
230
231 void
interpolate(double r,double s,const Ref<Vertex> & result,SCVector3 & dA,const Ref<Volume> & vol,double isoval)232 Triangle::interpolate(double r, double s,
233 const Ref<Vertex>&result, SCVector3&dA,
234 const Ref<Volume> &vol, double isoval)
235 {
236 // set up an initial dummy norm
237 SCVector3 norm(0.0);
238 result->set_normal(norm);
239
240 // initial guess
241 interpolate(r,s,result,dA);
242
243 // now refine that guess
244 SCVector3 trialpoint = result->point();
245 SCVector3 trialnorm = result->normal();
246 SCVector3 newpoint;
247 vol->solve(trialpoint,trialnorm,isoval,newpoint);
248 // compute the true normal
249 vol->set_x(newpoint);
250 if (vol->gradient_implemented()) {
251 vol->get_gradient(trialnorm);
252 }
253 trialnorm.normalize();
254 result->set_point(newpoint);
255 result->set_normal(trialnorm);
256 }
257
258 void
flip()259 Triangle::flip()
260 {
261 _orientation0 = _orientation0?0:1;
262 _orientation1 = _orientation1?0:1;
263 _orientation2 = _orientation2?0:1;
264 }
265
266 void
set_order(int order,const Ref<Volume> & vol,double isovalue)267 Triangle::set_order(int order, const Ref<Volume>&vol, double isovalue)
268 {
269 if (order > max_order) {
270 ExEnv::errn() << scprintf("Triangle::set_order: max_order = %d but order = %d\n",
271 max_order, order);
272 abort();
273 }
274
275 unsigned int i, j, k;
276
277 if (edge(0)->order() != order
278 ||edge(1)->order() != order
279 ||edge(2)->order() != order) {
280 ExEnv::errn() << "Triangle::set_order: edge order doesn't match" << endl;
281 abort();
282 }
283
284 _order = order;
285 delete[] _vertices;
286
287 _vertices = new Ref<Vertex>[TriInterpCoef::order_to_nvertex(_order)];
288
289 // fill in the corner vertices
290 _vertices[TriInterpCoef::ijk_to_index(_order, 0, 0)] = vertex(0);
291 _vertices[TriInterpCoef::ijk_to_index(0, 0, _order)] = vertex(1);
292 _vertices[TriInterpCoef::ijk_to_index(0, _order, 0)] = vertex(2);
293
294 // fill in the interior edge vertices
295 for (i = 1; i < _order; i++) {
296 j = _order - i;
297 _vertices[TriInterpCoef::ijk_to_index(0, i, j)]
298 = _edges[1]->interior_vertex(_orientation1?i:j);
299 _vertices[TriInterpCoef::ijk_to_index(j, 0, i)]
300 = _edges[0]->interior_vertex(_orientation0?i:j);
301 _vertices[TriInterpCoef::ijk_to_index(i, j, 0)]
302 = _edges[2]->interior_vertex(_orientation2?i:j);
303 }
304
305 const SCVector3& p0 = vertex(0)->point();
306 const SCVector3& p1 = vertex(1)->point();
307 const SCVector3& p2 = vertex(2)->point();
308 const SCVector3& norm0 = vertex(0)->normal();
309 const SCVector3& norm1 = vertex(1)->normal();
310 const SCVector3& norm2 = vertex(2)->normal();
311
312 for (i=0; i<=_order; i++) {
313 double I = (1.0*i)/_order;
314 for (j=0; j<=_order-i; j++) {
315 SCVector3 trialpoint;
316 SCVector3 trialnorm;
317 SCVector3 newpoint;
318 double J = (1.0*j)/_order;
319 k = _order - i - j;
320 if (!i || !j || !k) continue; // interior point check
321 double K = (1.0*k)/_order;
322 int index = TriInterpCoef::ijk_to_index(i,j,k);
323 // this get approximate vertices and normals
324 trialpoint = I*p0 + J*p1 + K*p2;
325 trialnorm = I*norm0 + J*norm1 + K*norm2;
326 // now refine that guess
327 vol->solve(trialpoint,trialnorm,isovalue,newpoint);
328 // compute the true normal
329 vol->set_x(newpoint);
330 if (vol->gradient_implemented()) {
331 vol->get_gradient(trialnorm);
332 }
333 trialnorm.normalize();
334 _vertices[index] = new Vertex(newpoint,trialnorm);
335 }
336 }
337 }
338
339 /////////////////////////////////////////////////////////////////////////
340 // TriangleIntegrator
341
342 static ClassDesc TriangleIntegrator_cd(
343 typeid(TriangleIntegrator),"TriangleIntegrator",1,"public DescribedClass",
344 0, create<TriangleIntegrator>, 0);
345
TriangleIntegrator(int order)346 TriangleIntegrator::TriangleIntegrator(int order):
347 _n(order)
348 {
349 _r = new double [_n];
350 _s = new double [_n];
351 _w = new double [_n];
352 coef_ = 0;
353 }
354
TriangleIntegrator(const Ref<KeyVal> & keyval)355 TriangleIntegrator::TriangleIntegrator(const Ref<KeyVal>& keyval)
356 {
357 _n = keyval->intvalue("n");
358 if (keyval->error() != KeyVal::OK) {
359 _n = 7;
360 }
361 _r = new double [_n];
362 _s = new double [_n];
363 _w = new double [_n];
364 coef_ = 0;
365 }
366
~TriangleIntegrator()367 TriangleIntegrator::~TriangleIntegrator()
368 {
369 delete[] _r;
370 delete[] _s;
371 delete[] _w;
372 clear_coef();
373 }
374
375 void
set_n(int n)376 TriangleIntegrator::set_n(int n)
377 {
378 delete[] _r;
379 delete[] _s;
380 delete[] _w;
381 _n = n;
382 _r = new double [_n];
383 _s = new double [_n];
384 _w = new double [_n];
385 clear_coef();
386 }
387
388 void
set_w(int i,double w)389 TriangleIntegrator::set_w(int i,double w)
390 {
391 _w[i] = w;
392 }
393
394 void
set_r(int i,double r)395 TriangleIntegrator::set_r(int i,double r)
396 {
397 _r[i] = r;
398 }
399
400 void
set_s(int i,double s)401 TriangleIntegrator::set_s(int i,double s)
402 {
403 _s[i] = s;
404 }
405
406 void
init_coef()407 TriangleIntegrator::init_coef()
408 {
409 int i, j;
410
411 clear_coef();
412
413 coef_ = new Ref<TriInterpCoef>*[Triangle::max_order];
414 for (i=0; i<Triangle::max_order; i++) {
415 coef_[i] = new Ref<TriInterpCoef>[_n];
416 for (j=0; j<_n; j++) {
417 TriInterpCoefKey key(i+1,_r[j],_s[j]);
418 coef_[i][j] = new TriInterpCoef(key);
419 }
420 }
421 }
422
423 void
clear_coef()424 TriangleIntegrator::clear_coef()
425 {
426 int i, j;
427
428 if (coef_) {
429 for (i=0; i<Triangle::max_order; i++) {
430 for (j=0; j<_n; j++) {
431 coef_[i][j] = 0;
432 }
433 delete[] coef_[i];
434 }
435 }
436 delete[] coef_;
437 coef_ = 0;
438 }
439
440 /////////////////////////////////////////////////////////////////////////
441 // GaussTriangleIntegrator
442
443 static ClassDesc GaussTriangleIntegrator_cd(
444 typeid(GaussTriangleIntegrator),"GaussTriangleIntegrator",1,"public TriangleIntegrator",
445 0, create<GaussTriangleIntegrator>, 0);
446
GaussTriangleIntegrator(const Ref<KeyVal> & keyval)447 GaussTriangleIntegrator::GaussTriangleIntegrator(const Ref<KeyVal>& keyval):
448 TriangleIntegrator(keyval)
449 {
450 ExEnv::out0() << "Created a GaussTriangleIntegrator with n = " << n() << endl;
451 init_rw(n());
452 init_coef();
453 }
454
GaussTriangleIntegrator(int order)455 GaussTriangleIntegrator::GaussTriangleIntegrator(int order):
456 TriangleIntegrator(order)
457 {
458 init_rw(n());
459 init_coef();
460 }
461
462 void
set_n(int n)463 GaussTriangleIntegrator::set_n(int n)
464 {
465 TriangleIntegrator::set_n(n);
466 init_rw(n);
467 init_coef();
468 }
469
470 void
init_rw(int order)471 GaussTriangleIntegrator::init_rw(int order)
472 {
473 if (order == 1) {
474 set_r(0, 1.0/3.0);
475 set_s(0, 1.0/3.0);
476 set_w(0, 1.0);
477 }
478 else if (order == 3) {
479 set_r(0, 1.0/6.0);
480 set_r(1, 2.0/3.0);
481 set_r(2, 1.0/6.0);
482 set_s(0, 1.0/6.0);
483 set_s(1, 1.0/6.0);
484 set_s(2, 2.0/3.0);
485 set_w(0, 1.0/3.0);
486 set_w(1, 1.0/3.0);
487 set_w(2, 1.0/3.0);
488 }
489 else if (order == 4) {
490 set_r(0, 1.0/3.0);
491 set_r(1, 1.0/5.0);
492 set_r(2, 3.0/5.0);
493 set_r(3, 1.0/5.0);
494 set_s(0, 1.0/3.0);
495 set_s(1, 1.0/5.0);
496 set_s(2, 1.0/5.0);
497 set_s(3, 3.0/5.0);
498 set_w(0, -27.0/48.0);
499 set_w(1, 25.0/48.0);
500 set_w(2, 25.0/48.0);
501 set_w(3, 25.0/48.0);
502 }
503 else if (order == 6) {
504 set_r(0, 0.091576213509771);
505 set_r(1, 0.816847572980459);
506 set_r(2, r(0));
507 set_r(3, 0.445948490915965);
508 set_r(4, 0.108103018168070);
509 set_r(5, r(3));
510 set_s(0, r(0));
511 set_s(1, r(0));
512 set_s(2, r(1));
513 set_s(3, r(3));
514 set_s(4, r(3));
515 set_s(5, r(4));
516 set_w(0, 0.109951743655322);
517 set_w(1, w(0));
518 set_w(2, w(0));
519 set_w(3, 0.223381589678011);
520 set_w(4, w(3));
521 set_w(5, w(3));
522 }
523 else if (order == 7) {
524 set_r(0, 1.0/3.0);
525 set_r(1, 0.101286507323456);
526 set_r(2, 0.797426985353087);
527 set_r(3, r(1));
528 set_r(4, 0.470142064105115);
529 set_r(5, 0.059715871789770);
530 set_r(6, r(4));
531 set_s(0, r(0));
532 set_s(1, r(1));
533 set_s(2, r(1));
534 set_s(3, r(2));
535 set_s(4, r(4));
536 set_s(5, r(4));
537 set_s(6, r(5));
538 set_w(0, 0.225);
539 set_w(1, 0.125939180544827);
540 set_w(2, w(1));
541 set_w(3, w(1));
542 set_w(4, 0.132394152788506);
543 set_w(5, w(4));
544 set_w(6, w(4));
545 }
546 else {
547 ExEnv::errn() << "GaussTriangleIntegrator: invalid order " << order << endl;
548 abort();
549 }
550
551 // scale the weights by the area of the unit triangle
552 for (int i=0; i<order; i++) {
553 set_w(i, w(i)*0.5);
554 }
555 }
556
~GaussTriangleIntegrator()557 GaussTriangleIntegrator::~GaussTriangleIntegrator()
558 {
559 }
560
561 /////////////////////////////////////////////////////////////////////////////
562
563 // Local Variables:
564 // mode: c++
565 // c-file-style: "CLJ"
566 // End:
567