source: issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp@ 5397

Last change on this file since 5397 was 5397, checked in by Mathieu Morlighem, 15 years ago

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