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

Last change on this file since 15061 was 15061, checked in by Eric.Larour, 12 years ago

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