source: issm/trunk-jpl/src/c/bamg/BamgVertex.cpp@ 18064

Last change on this file since 18064 was 18064, checked in by Mathieu Morlighem, 11 years ago

DEL: removed all FUNCTIONs
File size: 6.4 KB
Line 
1#include <cstdio>
2#include <cstring>
3#include <cmath>
4#include <ctime>
5
6#include "./bamgobjects.h"
7#include "../shared/shared.h"
8#include "./det.h"
9
10namespace bamg {
11
12 /*Methods*/
13 void BamgVertex::Echo(void){/*{{{*/
14
15 _printf_("Vertex:\n");
16 _printf_(" integer coordinates i.x: " << i.x << ", i.y: " << i.y << "\n");
17 _printf_(" Euclidean coordinates r.x: " << r.x << ", r.y: " << r.y << "\n");
18 _printf_(" ReferenceNumber = " << ReferenceNumber << "\n");
19 m.Echo();
20
21 return;
22 }
23 /*}}}*/
24 int BamgVertex::GetReferenceNumber() const { /*{{{*/
25 return ReferenceNumber;
26 }
27 /*}}}*/
28 void BamgVertex::MetricFromHessian(const double Hxx,const double Hyx, const double Hyy,const double smin,const double smax,const double s,double err,BamgOpts* bamgopts){/*{{{*/
29 /*Compute Metric from Hessian*/
30
31 /*get options*/
32 double power=(bamgopts->power);
33 double anisomax=(bamgopts->anisomax);
34 double CutOff=bamgopts->cutoff;
35 double hmin=(bamgopts->hmin);
36 double hmax=(bamgopts->hmax);
37 double coef=bamgopts->coeff;
38 int Metrictype=(bamgopts->Metrictype);
39
40 /*Intermediary*/
41 double ci;
42
43 /*compute multiplicative coefficient depending on Metric Type (2/9 because it is 2d)*/
44
45 //Absolute Error
46 /*
47 * 2 1
48 *Metric M = --- ------------ Abs(Hessian)
49 * 9 err * coeff^2
50 */
51 if (Metrictype==0){
52 ci= 2.0/9.0 * 1/(err*coef*coef);
53 }
54
55 //Relative Error
56 /*
57 * 2 1 Abs(Hessian)
58 *Metric M = --- ------------ ----------------------
59 * 9 err * coeff^2 max( |s| , cutoff*max(|s|) )
60 *
61 */
62 else if (Metrictype==1){
63 ci= 2.0/9.0 * 1/(err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
64 }
65
66 //Rescaled absolute error
67 /*
68 * 2 1 Abs(Hessian)
69 *Metric M = --- ------------ ----------------------
70 * 9 err * coeff^2 (smax-smin)
71 */
72 else if (Metrictype==2){
73 ci= 2.0/9.0 * 1/(err*coef*coef) * 1/(smax-smin);
74 }
75 else{
76 _error_("Metrictype " << Metrictype << " not supported yet (use 0,1 or 2(default))");
77 }
78
79 //initialize metric Miv with ci*H
80 Metric Miv(Hxx*ci,Hyx*ci,Hyy*ci);
81
82 //Get eigen values and vectors of Miv
83 EigenMetric Vp(Miv);
84
85 //move eigen valuse to their absolute values
86 Vp.Abs();
87
88 //Apply a power if requested by user
89 if(power!=1.0) Vp.pow(power);
90
91 //modify eigen values according to hmin and hmax
92 Vp.Maxh(hmax);
93 Vp.Minh(hmin);
94
95 //Bound anisotropy by 1/(anisomax)^2
96 Vp.BoundAniso2(1/(anisomax*anisomax));
97
98 //rebuild Metric from Vp
99 Metric MVp(Vp);
100
101 //Apply Metric to vertex
102 m.IntersectWith(MVp);
103
104 }
105 /*}}}*/
106 long BamgVertex::Optim(int i,int koption){ /*{{{*/
107 long ret=0;
108 if ( t && (IndexInTriangle >= 0 ) && (IndexInTriangle <3) ){
109 ret = t->Optim(IndexInTriangle,koption);
110 if(!i){
111 t =0; // for no future optime
112 IndexInTriangle= 0;
113 }
114 }
115 return ret;
116 }
117 /*}}}*/
118 double BamgVertex::Smoothing(Mesh &Th,const Mesh &BTh,Triangle* &tstart ,double omega){/*{{{*/
119 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/
120
121 BamgVertex* s=this;
122 BamgVertex &vP = *s,vPsave=vP;
123
124 Triangle* tbegin= t , *tria = t , *ttc;
125
126 int k=0,kk=0,j = EdgesVertexTriangle[IndexInTriangle][0],jc;
127 R2 P(s->r),PNew(0,0);
128 do {
129 k++;
130
131 if (!tria->Hidden(j)){
132 BamgVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
133
134 R2 Q = vQ,QP(P-Q);
135 double lQP = LengthInterpole(vP,vQ,QP);
136 PNew += Q+QP/Max(lQP,1e-20);
137 kk ++;
138 }
139 ttc = tria->TriangleAdj(j);
140 jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
141 tria = ttc;
142 j = NextEdge[jc];
143 if (k>=2000){
144 _error_("k>=2000 (Maximum number of iterations reached)");
145 }
146 } while ( tbegin != tria);
147 if (kk<4) return 0;
148 PNew = PNew/(double)kk;
149 R2 Xmove((PNew-P)*omega);
150 PNew = P+Xmove;
151 double delta=Norme2_2(Xmove);
152
153 Icoor2 deta[3];
154 I2 IBTh = BTh.R2ToI2(PNew);
155
156 tstart=BTh.TriangleFindFromCoord(IBTh,deta,tstart);
157
158 if (tstart->det <0){ // outside
159 double ba,bb;
160 AdjacentTriangle edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ;
161 tstart = edge;
162 vP.m= Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));
163 }
164 else { // inside
165 double aa[3];
166 double s = deta[0]+deta[1]+deta[2];
167 aa[0]=deta[0]/s;
168 aa[1]=deta[1]/s;
169 aa[2]=deta[2]/s;
170 vP.m = Metric(aa,(*tstart)[0],(*tstart)[1],(*tstart)[2]);
171 }
172
173 // recompute the det of the triangle
174 vP.r = PNew;
175
176 vP.i = Th.R2ToI2(PNew);
177
178 BamgVertex vPnew = vP;
179
180 int ok=1;
181 int loop=1;
182 k=0;
183 while (ok){
184 ok =0;
185 do {
186 k++;
187 double detold = tria->det;
188 tria->det = bamg::det( (*tria)[0],(*tria)[1] ,(*tria)[2]);
189 if (loop) {
190 BamgVertex *v0,*v1,*v2,*v3;
191 if (tria->det<0) ok =1;
192 else if (tria->Quadrangle(v0,v1,v2,v3))
193 {
194 vP = vPsave;
195 double qold =QuadQuality(*v0,*v1,*v2,*v3);
196 vP = vPnew;
197 double qnew =QuadQuality(*v0,*v1,*v2,*v3);
198 if (qnew<qold) ok = 1;
199 }
200 else if ( (double)tria->det < detold/2 ) ok=1;
201
202 }
203 tria->SetUnMarkUnSwap(0);
204 tria->SetUnMarkUnSwap(1);
205 tria->SetUnMarkUnSwap(2);
206 ttc = tria->TriangleAdj(j);
207 jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
208 tria = ttc;
209 j = NextEdge[jc];
210 if (k>=2000){
211 _error_("k>=2000");
212 }
213 }while ( tbegin != tria);
214
215 if (ok && loop) vP=vPsave; // no move
216 loop=0;
217 }
218 return delta;
219 }
220 /*}}}*/
221
222 /*Intermediary*/
223 double QuadQuality(const BamgVertex & a,const BamgVertex &b,const BamgVertex &c,const BamgVertex &d) {/*{{{*/
224 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshQuad.cpp/QuadQuality)*/
225
226 // calcul de 4 angles --
227 R2 A((R2)a),B((R2)b),C((R2)c),D((R2)d);
228 R2 AB(B-A),BC(C-B),CD(D-C),DA(A-D);
229 // Move(A),Line(B),Line(C),Line(D),Line(A);
230 const Metric & Ma = a;
231 const Metric & Mb = b;
232 const Metric & Mc = c;
233 const Metric & Md = d;
234
235 double lAB=Norme2(AB);
236 double lBC=Norme2(BC);
237 double lCD=Norme2(CD);
238 double lDA=Norme2(DA);
239 AB /= lAB; BC /= lBC; CD /= lCD; DA /= lDA;
240 // version aniso
241 double cosDAB= Ma(DA,AB)/(Ma(DA)*Ma(AB)),sinDAB= Det(DA,AB);
242 double cosABC= Mb(AB,BC)/(Mb(AB)*Mb(BC)),sinABC= Det(AB,BC);
243 double cosBCD= Mc(BC,CD)/(Mc(BC)*Mc(CD)),sinBCD= Det(BC,CD);
244 double cosCDA= Md(CD,DA)/(Md(CD)*Md(DA)),sinCDA= Det(CD,DA);
245 double sinmin=Min(Min(sinDAB,sinABC),Min(sinBCD,sinCDA));
246 if (sinmin<=0) return sinmin;
247 return 1.0-Max(Max(Abs(cosDAB),Abs(cosABC)),Max(Abs(cosBCD),Abs(cosCDA)));
248 }
249 /*}}}*/
250
251}
Note: See TracBrowser for help on using the repository browser.