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

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

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