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

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

moved Number to GetId

File size: 3.1 KB
RevLine 
[3913]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::R1tg{{{1*/
19 double GeometricalEdge::R1tg(double theta,R2 & t) const{
20 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/R1tg)*/
21 // 1/R of radius of cuvature
22
23 R2 A=v[0]->r,B=v[1]->r;
24 double dca,dcb,dcta,dctb;
25 double ddca,ddcb,ddcta,ddctb;
26 double tt = theta*theta;
27
28 //check theta
29 ISSMASSERT(theta>=0 && theta<=1);
30
31 if (TgA()){
32 if (TgB()){
33 // Tangent A and B provided:
34 // interpolation d'hermite
35 dcb = 6*theta*(1-theta);
36 ddcb = 6*(1-2*theta);
37 dca = -dcb;
38 ddca = -ddcb;
39 dcta = (3*theta - 4)*theta + 1;
40 ddcta=6*theta-4;
41 dctb = 3*tt - 2*theta;
42 ddctb = 6*theta-2;
43 }
44 else {
45 //Tangent A provided but tangent B not provided
46 // 1-t*t, t-t*t, t*t
47 double t = theta;
48 dcb = 2*t;
49 ddcb = 2;
50 dca = -dcb;
51 ddca = -2;
52 dcta = 1-dcb;
53 ddcta = -ddcb;
54 dctb=0;
55 ddctb=0;
56 }
57 }
58 else{
59 if (TgB()){
60 //Tangent B provided but tangent A not provided
61 double t = 1-theta;
62 dca = -2*t;
63 ddca = 2;
64 dcb = -dca;
65 ddcb = -2;
66 dctb = 1+dca;
67 ddctb= ddca;
68 dcta =0;
69 ddcta =0;
70 }
71 else {
72 //Neither thangent A nor tangent B provided
73 // lagrange P1
74 t=B-A;
75 return 0;
76 }
77 }
78 R2 d = A*dca + B*dcb + tg[0]* dcta + tg[1] * dctb;
79 R2 dd = A*ddca + B*ddcb + tg[0]* ddcta + tg[1] * ddctb;
80 double d2=(d,d);
81 double sd2 = sqrt(d2);
82 t=d;
83 if(d2>1.0e-20){
84 t/=sd2;
85 return Abs(Det(d,dd))/(d2*sd2);
86 }
87 else return 0;
88 }
89 /*}}}1*/
90 /*FUNCTION GeometricalEdge::F{{{1*/
91 R2 GeometricalEdge::F(double theta) const{
92 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/F)*/
93 // parametrization of the curve edge
94
95 R2 A=v[0]->r,B=v[1]->r;
96 double ca,cb,cta,ctb;
97 if ( theta<-1e-12){
98 ISSMERROR("theta<-1e-12");
99 }
100 if ( theta>1+1e-12){
101 ISSMERROR("theta>1+1e-12");
102 }
103 if (TgA())
104 if (TgB()) // interpolation d'hermite
105 { cb = theta*theta*(3-2*theta);
106 ca = 1-cb;
107 cta = (1-theta)*(1-theta)*theta;
108 ctb = (theta-1)*theta*theta ;
109 }
110 else { // 1-t*t, t-t*t, t*t
111 double t = theta;
112 cb = t*t;
113 ca = 1-cb;
114 cta= t-cb;
115 ctb=0;
116 }
117 else
118 if (TgB()){
119 double t = 1-theta;
120 ca = t*t;
121 cb = 1-ca;
122 ctb= -t+ca;
123 cta=0;
124 }
125 else {
126 ca =(1-theta),cb = theta,cta=ctb=0; // lagrange P1
127 }
128 return A*ca + B*cb + tg[0]* cta + tg[1] * ctb;
129 }
130 /*}}}1*/
131 /*FUNCTION GeometricalEdge::Set {{{1*/
132 void GeometricalEdge::Set(const GeometricalEdge & rec,const Geometry & Gh ,Geometry & GhNew){
133 *this = rec;
[5149]134 v[0] = GhNew.vertices + Gh.GetId(v[0]);
135 v[1] = GhNew.vertices + Gh.GetId(v[1]);
136 if (Adj[0]) Adj[0] = GhNew.edges + Gh.GetId(Adj[0]);
137 if (Adj[1]) Adj[1] = GhNew.edges + Gh.GetId(Adj[1]);
[3913]138 }
139 /*}}}*/
140
141}
Note: See TracBrowser for help on using the repository browser.