source: ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.cc@ 91c409

Candidate_v1.7.0 stable
Last change on this file since 91c409 was 860145, checked in by Frederik Heber <heber@…>, 9 years ago

Merge commit '0b990dfaa8c6007a996d030163a25f7f5fc8a7e7' as 'ThirdParty/mpqc_open'
  • Property mode set to 100644
File size: 15.4 KB
RevLine 
[0b990d]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
37using namespace std;
38using 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
71Triangle::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
133Triangle::~Triangle()
134{
135 if (_vertices) delete[] _vertices;
136}
137
138Ref<Vertex>
139Triangle::vertex(int i)
140{
141 return _edges[i]->vertex(orientation(i));
142}
143
144unsigned int
145Triangle::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
152int
153Triangle::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
161double 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
173void 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
178void Triangle::add_edges(std::set<Ref<Edge> >&set)
179{
180 for (int i=0; i<3; i++) set.insert(_edges[i]);
181}
182
183void
184Triangle::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
191void
192Triangle::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
231void
232Triangle::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
258void
259Triangle::flip()
260{
261 _orientation0 = _orientation0?0:1;
262 _orientation1 = _orientation1?0:1;
263 _orientation2 = _orientation2?0:1;
264}
265
266void
267Triangle::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
342static ClassDesc TriangleIntegrator_cd(
343 typeid(TriangleIntegrator),"TriangleIntegrator",1,"public DescribedClass",
344 0, create<TriangleIntegrator>, 0);
345
346TriangleIntegrator::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
355TriangleIntegrator::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
367TriangleIntegrator::~TriangleIntegrator()
368{
369 delete[] _r;
370 delete[] _s;
371 delete[] _w;
372 clear_coef();
373}
374
375void
376TriangleIntegrator::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
388void
389TriangleIntegrator::set_w(int i,double w)
390{
391 _w[i] = w;
392}
393
394void
395TriangleIntegrator::set_r(int i,double r)
396{
397 _r[i] = r;
398}
399
400void
401TriangleIntegrator::set_s(int i,double s)
402{
403 _s[i] = s;
404}
405
406void
407TriangleIntegrator::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
423void
424TriangleIntegrator::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
443static ClassDesc GaussTriangleIntegrator_cd(
444 typeid(GaussTriangleIntegrator),"GaussTriangleIntegrator",1,"public TriangleIntegrator",
445 0, create<GaussTriangleIntegrator>, 0);
446
447GaussTriangleIntegrator::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
455GaussTriangleIntegrator::GaussTriangleIntegrator(int order):
456 TriangleIntegrator(order)
457{
458 init_rw(n());
459 init_coef();
460}
461
462void
463GaussTriangleIntegrator::set_n(int n)
464{
465 TriangleIntegrator::set_n(n);
466 init_rw(n);
467 init_coef();
468}
469
470void
471GaussTriangleIntegrator::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
557GaussTriangleIntegrator::~GaussTriangleIntegrator()
558{
559}
560
561/////////////////////////////////////////////////////////////////////////////
562
563// Local Variables:
564// mode: c++
565// c-file-style: "CLJ"
566// End:
Note: See TracBrowser for help on using the repository browser.