source: issm/trunk/src/c/objects/Bamg/GeomEdge.cpp@ 6412

Last change on this file since 6412 was 6412, checked in by Mathieu Morlighem, 14 years ago

moved ISSMERROR to _error_, ISSMASSERT to _assert_ and ISSMPRINTF to _printf_
File size: 4.3 KB
Line 
1#include <cstdio>
2#include <string.h>
3#include <cmath>
4
5#include "../../include/include.h"
6#include "../../shared/Exceptions/exceptions.h"
7
8#include "GeomEdge.h"
9#include "Geometry.h"
10
11using namespace std;
12
13namespace bamg {
14
15 /*Constructor/Destructor*/
16
17 /*Methods*/
18 /*FUNCTION GeomEdge::Cracked{{{1*/
19 int GeomEdge::Cracked() const {
20 return type &1;
21 }/*}}}*/
22 /*FUNCTION GeomEdge::F{{{1*/
23 R2 GeomEdge::F(double theta) const{
24 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/F)*/
25 // parametrization of the curve edge
26
27 R2 A=v[0]->r,B=v[1]->r;
28 double ca,cb,cta,ctb;
29
30 //Check that theta is in [0 1]
31 _assert_(theta>-1e-12 && theta<1+1e-12);
32
33 if (TgA()){
34 if (TgB()){ //Hermite interpolation
35 cb = theta*theta*(3-2*theta);
36 ca = 1-cb;
37 cta = (1-theta)*(1-theta)*theta;
38 ctb = (theta-1)*theta*theta ;
39 }
40 else {
41 double t = theta;
42 cb = t*t;
43 ca = 1-cb;
44 cta= t-cb;
45 ctb=0;
46 }
47 }
48 else{
49 if (TgB()){
50 double t = 1-theta;
51 ca = t*t;
52 cb = 1-ca;
53 ctb= -t+ca;
54 cta=0;
55 }
56 else { // lagrange P1
57 ca =(1-theta);
58 cb = theta;
59 cta=ctb=0;
60 }
61 }
62 return A*ca + B*cb + tg[0]*cta + tg[1]*ctb;
63 }
64 /*}}}1*/
65 /*FUNCTION GeomEdge::Mark{{{1*/
66 int GeomEdge::Mark() const {
67 return type &16;
68 }/*}}}*/
69 /*FUNCTION GeomEdge::R1tg{{{1*/
70 double GeomEdge::R1tg(double theta,R2 & t) const{
71 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/R1tg)*/
72 // 1/R of radius of cuvature
73
74 R2 A=v[0]->r,B=v[1]->r;
75 double dca,dcb,dcta,dctb;
76 double ddca,ddcb,ddcta,ddctb;
77 double tt = theta*theta;
78
79 //check theta
80 _assert_(theta>=0 && theta<=1);
81
82 if (TgA()){
83 if (TgB()){
84 // Tangent A and B provided:
85 // interpolation d'hermite
86 dcb = 6*theta*(1-theta);
87 ddcb = 6*(1-2*theta);
88 dca = -dcb;
89 ddca = -ddcb;
90 dcta = (3*theta - 4)*theta + 1;
91 ddcta=6*theta-4;
92 dctb = 3*tt - 2*theta;
93 ddctb = 6*theta-2;
94 }
95 else {
96 //Tangent A provided but tangent B not provided
97 // 1-t*t, t-t*t, t*t
98 double t = theta;
99 dcb = 2*t;
100 ddcb = 2;
101 dca = -dcb;
102 ddca = -2;
103 dcta = 1-dcb;
104 ddcta = -ddcb;
105 dctb=0;
106 ddctb=0;
107 }
108 }
109 else{
110 if (TgB()){
111 //Tangent B provided but tangent A not provided
112 double t = 1-theta;
113 dca = -2*t;
114 ddca = 2;
115 dcb = -dca;
116 ddcb = -2;
117 dctb = 1+dca;
118 ddctb= ddca;
119 dcta =0;
120 ddcta =0;
121 }
122 else {
123 //Neither thangent A nor tangent B provided
124 // lagrange P1
125 t=B-A;
126 return 0;
127 }
128 }
129 R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb;
130 R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;
131 double d2=(d,d);
132 double sd2 = sqrt(d2);
133 t=d;
134 if(d2>1.0e-20){
135 t/=sd2;
136 return Abs(Det(d,dd))/(d2*sd2);
137 }
138 else return 0;
139 }
140 /*}}}1*/
141 /*FUNCTION GeomEdge::Required{{{1*/
142 int GeomEdge::Required() {
143 return type &64;
144 }/*}}}*/
145 /*FUNCTION GeomEdge::Set {{{1*/
146 void GeomEdge::Set(const GeomEdge & rec,const Geometry & Gh ,Geometry & GhNew){
147 *this = rec;
148 v[0] = GhNew.vertices + Gh.GetId(v[0]);
149 v[1] = GhNew.vertices + Gh.GetId(v[1]);
150 if (Adj[0]) Adj[0] = GhNew.edges + Gh.GetId(Adj[0]);
151 if (Adj[1]) Adj[1] = GhNew.edges + Gh.GetId(Adj[1]);
152 }
153 /*}}}*/
154 /*FUNCTION GeomEdge::SetCracked{{{1*/
155 void GeomEdge::SetCracked() {
156 type |= 1;/*=>1st digit to 1*/
157 }/*}}}*/
158 /*FUNCTION GeomEdge::SetTgA{{{1*/
159 void GeomEdge::SetTgA() {
160 type |=4; /*=>2d digit to 1*/
161 }/*}}}*/
162 /*FUNCTION GeomEdge::SetTgB{{{1*/
163 void GeomEdge::SetTgB() {
164 type |=8; /*=> 3d digit to 1*/
165 }/*}}}*/
166 /*FUNCTION GeomEdge::SetMark{{{1*/
167 void GeomEdge::SetMark() {
168 type |=16;/*=> 4th digiy to 1*/
169 }/*}}}*/
170 /*FUNCTION GeomEdge::SetUnMark{{{1*/
171 void GeomEdge::SetUnMark() {
172 type &= 1007 /* 1023-16 = 000111110111 => 4th digit to 0*/;
173 }/*}}}*/
174 /*FUNCTION GeomEdge::SetRequired{{{1*/
175 void GeomEdge::SetRequired() {
176 type |= 64;/*=>6th digit to 1*/
177 }/*}}}*/
178 /*FUNCTION GeomEdge::Tg{{{1*/
179 int GeomEdge::Tg(int i) const {
180 return i==0 ? TgA() : TgB();
181 }/*}}}*/
182 /*FUNCTION GeomEdge::TgA{{{1*/
183 int GeomEdge::TgA() const {
184 return type &4;
185 }/*}}}*/
186 /*FUNCTION GeomEdge::TgB{{{1*/
187 int GeomEdge::TgB() const {
188 return type &8;
189 }/*}}}*/
190}
Note: See TracBrowser for help on using the repository browser.