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

Last change on this file was 22253, checked in by tsantos, 7 years ago

CHG: adding verbose in src/c/bamg/Mesh.cpp
File size: 5.2 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 /*Constructor/Destructor*/
13 BamgVertex::BamgVertex(){ /*{{{*/
14 this->PreviousNumber = 0;
15 }/*}}}*/
16
17 /*Methods*/
18 void BamgVertex::Echo(void){/*{{{*/
19
20 _printf_("Vertex:\n");
21 _printf_(" integer coordinates i.x: " << i.x << ", i.y: " << i.y << "\n");
22 _printf_(" Euclidean coordinates r.x: " << r.x << ", r.y: " << r.y << "\n");
23 _printf_(" ReferenceNumber = " << ReferenceNumber << "\n");
24 _printf_(" PreviousNumber = " << PreviousNumber << "\n");
25 m.Echo();
26
27 return;
28 }
29 /*}}}*/
30 int BamgVertex::GetReferenceNumber() const { /*{{{*/
31 return ReferenceNumber;
32 }
33 /*}}}*/
34 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){/*{{{*/
35 /*Compute Metric from Hessian*/
36
37 /*get options*/
38 double power=(bamgopts->power);
39 double anisomax=(bamgopts->anisomax);
40 double CutOff=bamgopts->cutoff;
41 double hmin=(bamgopts->hmin);
42 double hmax=(bamgopts->hmax);
43 double coef=bamgopts->coeff;
44 int Metrictype=(bamgopts->Metrictype);
45
46 /*Intermediary*/
47 double ci;
48
49 /*compute multiplicative coefficient depending on Metric Type (2/9 because it is 2d)*/
50
51 //Absolute Error
52 /*
53 * 2 1
54 *Metric M = --- ------------ Abs(Hessian)
55 * 9 err * coeff^2
56 */
57 if (Metrictype==0){
58 ci= 2.0/9.0 * 1/(err*coef*coef);
59 }
60
61 //Relative Error
62 /*
63 * 2 1 Abs(Hessian)
64 *Metric M = --- ------------ ----------------------
65 * 9 err * coeff^2 max( |s| , cutoff*max(|s|) )
66 *
67 */
68 else if (Metrictype==1){
69 ci= 2.0/9.0 * 1/(err*coef*coef) * 1/Max( Abs(s), CutOff*(Max(Abs(smin),Abs(smax))));
70 }
71
72 //Rescaled absolute error
73 /*
74 * 2 1 Abs(Hessian)
75 *Metric M = --- ------------ ----------------------
76 * 9 err * coeff^2 (smax-smin)
77 */
78 else if (Metrictype==2){
79 ci= 2.0/9.0 * 1/(err*coef*coef) * 1/(smax-smin);
80 }
81 else{
82 _error_("Metrictype " << Metrictype << " not supported yet (use 0,1 or 2(default))");
83 }
84
85 //initialize metric Miv with ci*H
86 Metric Miv(Hxx*ci,Hyx*ci,Hyy*ci);
87
88 //Get eigen values and vectors of Miv
89 EigenMetric Vp(Miv);
90
91 //move eigen valuse to their absolute values
92 Vp.Abs();
93
94 //Apply a power if requested by user
95 if(power!=1.0) Vp.pow(power);
96
97 //modify eigen values according to hmin and hmax
98 Vp.Maxh(hmax);
99 Vp.Minh(hmin);
100
101 //Bound anisotropy by 1/(anisomax)^2
102 Vp.BoundAniso2(1/(anisomax*anisomax));
103
104 //rebuild Metric from Vp
105 Metric MVp(Vp);
106
107 //Apply Metric to vertex
108 m.IntersectWith(MVp);
109
110 }
111 /*}}}*/
112 long BamgVertex::Optim(int i,int koption){ /*{{{*/
113 long ret=0;
114 if ( t && (IndexInTriangle >= 0 ) && (IndexInTriangle <3) ){
115 ret = t->Optim(IndexInTriangle,koption);
116 if(i==0){
117 t =0; // for no future optim
118 IndexInTriangle= 0;
119 }
120 }
121 return ret;
122 }
123 /*}}}*/
124 double BamgVertex::Smoothing(Mesh &Th,Mesh &BTh,Triangle* &tstart ,double omega){/*{{{*/
125 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/Smoothing)*/
126
127 BamgVertex* s=this;
128 BamgVertex &vP = *s,vPsave=vP;
129
130 Triangle* tbegin= t , *tria = t , *ttc;
131
132 int k=0,kk=0,j = EdgesVertexTriangle[IndexInTriangle][0],jc;
133 R2 P(s->r),PNew(0,0);
134 do {
135 k++;
136
137 if (!tria->Hidden(j)){
138 BamgVertex &vQ = (*tria)[VerticesOfTriangularEdge[j][0]];
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){
150 _error_("k>=2000 (Maximum number of iterations reached)");
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 long long deta[3];
160 I2 IBTh = BTh.R2ToI2(PNew);
161
162 tstart=BTh.TriangleFindFromCoord(IBTh,deta,tstart);
163
164 if (tstart->det <0){ // outside
165 double ba,bb;
166 AdjacentTriangle edge= CloseBoundaryEdge(IBTh,tstart,ba,bb) ;
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
182 vP.i = Th.R2ToI2(PNew);
183
184 BamgVertex vPnew = vP;
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) {
196 if (tria->det<0) ok =1;
197 else if ( (double)tria->det < detold/2 ) ok=1;
198 }
199 tria->SetUnMarkUnSwap(0);
200 tria->SetUnMarkUnSwap(1);
201 tria->SetUnMarkUnSwap(2);
202 ttc = tria->TriangleAdj(j);
203 jc = NextEdge[tria->NuEdgeTriangleAdj(j)];
204 tria = ttc;
205 j = NextEdge[jc];
206 if (k>=2000){
207 _error_("k>=2000");
208 }
209 }while ( tbegin != tria);
210
211 if (ok && loop) vP=vPsave; // no move
212 loop=0;
213 }
214 return delta;
215 }
216 /*}}}*/
217}
Note: See TracBrowser for help on using the repository browser.