[3913] | 1 | #include <cstdio>
|
---|
| 2 | #include <cstring>
|
---|
| 3 | #include <cmath>
|
---|
| 4 | #include <ctime>
|
---|
| 5 |
|
---|
| 6 | #include "../objects.h"
|
---|
| 7 | #include "../../include/include.h"
|
---|
| 8 | #include "../../shared/Exceptions/exceptions.h"
|
---|
| 9 |
|
---|
| 10 | namespace bamg {
|
---|
| 11 |
|
---|
| 12 | static const Direction NoDirOfSearch=Direction();
|
---|
| 13 |
|
---|
| 14 | /*Constructors/Destructors*/
|
---|
[12365] | 15 | /*FUNCTION Geometry::Geometry(){{{*/
|
---|
[5120] | 16 | Geometry::Geometry(){
|
---|
| 17 | Init();
|
---|
| 18 | }
|
---|
[5130] | 19 | /*}}}*/
|
---|
[12365] | 20 | /*FUNCTION Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){{{*/
|
---|
[5120] | 21 | Geometry::Geometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
|
---|
| 22 | Init();
|
---|
| 23 | ReadGeometry(bamggeom,bamgopts);
|
---|
[5397] | 24 | PostRead();
|
---|
[5120] | 25 | }
|
---|
[5130] | 26 | /*}}}*/
|
---|
[12365] | 27 | /*FUNCTION Geometry::Geometry(const Geometry & Gh) (COPY operator){{{*/
|
---|
[3913] | 28 | Geometry::Geometry(const Geometry & Gh) {
|
---|
| 29 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Geometry)*/
|
---|
| 30 |
|
---|
| 31 | long i;
|
---|
| 32 | *this = Gh;
|
---|
| 33 | NbRef =0;
|
---|
| 34 | quadtree=0;
|
---|
[5573] | 35 | vertices = nbv ? new GeomVertex[nbv] : NULL;
|
---|
| 36 | edges = nbe ? new GeomEdge[nbe]:NULL;
|
---|
[5148] | 37 | curves= nbcurves ? new Curve[nbcurves]:NULL;
|
---|
[5573] | 38 | subdomains = nbsubdomains ? new GeomSubDomain[nbsubdomains]:NULL;
|
---|
[3913] | 39 | for (i=0;i<nbe;i++)
|
---|
| 40 | edges[i].Set(Gh.edges[i],Gh,*this);
|
---|
[5148] | 41 | for (i=0;i<nbcurves;i++)
|
---|
[3913] | 42 | curves[i].Set(Gh.curves[i],Gh,*this);
|
---|
[5148] | 43 | for (i=0;i<nbsubdomains;i++)
|
---|
[3913] | 44 | subdomains[i].Set(Gh.subdomains[i],Gh,*this);
|
---|
| 45 | }
|
---|
[12365] | 46 | /*}}}*/
|
---|
| 47 | /*FUNCTION Geometry::~Geometry(){{{*/
|
---|
[3913] | 48 | Geometry::~Geometry() {
|
---|
| 49 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/~Geometry)*/
|
---|
[12509] | 50 | if(NbRef>0){ _printString_("Trying to delete geometry and NbRef>0, probably due to an error"); return;}
|
---|
[5208] | 51 | if(vertices) delete [] vertices; vertices=0;
|
---|
| 52 | if(edges) delete [] edges; edges=0;
|
---|
| 53 | if(quadtree) delete quadtree; quadtree=0;
|
---|
| 54 | if(curves) delete [] curves; curves=0;nbcurves=0;
|
---|
[3913] | 55 | if(subdomains) delete [] subdomains;subdomains=0;
|
---|
[5120] | 56 | Init();
|
---|
[3913] | 57 | }
|
---|
[12365] | 58 | /*}}}*/
|
---|
[3913] | 59 |
|
---|
| 60 | /*IO*/
|
---|
[12365] | 61 | /*FUNCTION Geometry::ReadGeometry{{{*/
|
---|
[3913] | 62 | void Geometry::ReadGeometry(BamgGeom* bamggeom,BamgOpts* bamgopts){
|
---|
| 63 |
|
---|
| 64 | int verbose;
|
---|
[5120] | 65 | nbv=0;
|
---|
| 66 | nbe=0;
|
---|
[5148] | 67 | nbcurves=0;
|
---|
[3913] | 68 |
|
---|
| 69 | double Hmin = HUGE_VAL;// the infinie value
|
---|
[5143] | 70 | int i,j,k,n,i0,i1,i2,i3;
|
---|
[3913] | 71 |
|
---|
[5120] | 72 | /*initialize some variables*/
|
---|
[5143] | 73 | verbose= bamgopts->verbose;
|
---|
| 74 | nbv = bamggeom->VerticesSize[0];
|
---|
| 75 | nbe = bamggeom->EdgesSize[0];
|
---|
[3913] | 76 |
|
---|
| 77 | //some checks
|
---|
[13036] | 78 | if (bamggeom->Vertices==NULL) _error_("the domain provided does not contain any vertex");
|
---|
| 79 | if (bamggeom->Edges==NULL) _error_("the domain provided does not contain any edge");
|
---|
[3913] | 80 |
|
---|
| 81 | //Vertices
|
---|
| 82 | if (bamggeom->Vertices){
|
---|
[12509] | 83 | if(verbose>5) _printLine_(" processing Vertices");
|
---|
[13036] | 84 | if (bamggeom->VerticesSize[1]!=3) _error_("Vertices should have 3 columns");
|
---|
[5573] | 85 | vertices = new GeomVertex[nbv];
|
---|
[3913] | 86 | for (i=0;i<nbv;i++) {
|
---|
| 87 | vertices[i].r.x=(double)bamggeom->Vertices[i*3+0];
|
---|
| 88 | vertices[i].r.y=(double)bamggeom->Vertices[i*3+1];
|
---|
| 89 | vertices[i].ReferenceNumber=(long)bamggeom->Vertices[i*3+2];
|
---|
| 90 | vertices[i].DirOfSearch=NoDirOfSearch;
|
---|
| 91 | vertices[i].color =0;
|
---|
[5124] | 92 | vertices[i].type=0;
|
---|
[3913] | 93 | }
|
---|
[5120] | 94 | /*find domain extrema (pmin,pmax) that will define the square box used for by the quadtree*/
|
---|
[3913] | 95 | pmin = vertices[0].r;
|
---|
| 96 | pmax = vertices[0].r;
|
---|
| 97 | for (i=0;i<nbv;i++) {
|
---|
| 98 | pmin.x = Min(pmin.x,vertices[i].r.x);
|
---|
| 99 | pmin.y = Min(pmin.y,vertices[i].r.y);
|
---|
| 100 | pmax.x = Max(pmax.x,vertices[i].r.x);
|
---|
| 101 | pmax.y = Max(pmax.y,vertices[i].r.y);
|
---|
| 102 | }
|
---|
[5120] | 103 | /*Offset pmin and pmax to avoid round-off errors*/
|
---|
| 104 | R2 offset = (pmax-pmin)*0.05;
|
---|
| 105 | pmin -= offset;
|
---|
| 106 | pmax += offset;
|
---|
| 107 | /*coefIcoor is the coefficient used for integer coordinates:
|
---|
| 108 | * (x-pmin.x)
|
---|
| 109 | * Icoor x = (2^30 -1) ------------
|
---|
| 110 | * D
|
---|
| 111 | * where D is the longest side of the domain (direction x or y)
|
---|
| 112 | * so that (x-pmin.x)/D is in ]0 1[
|
---|
[5177] | 113 | *
|
---|
| 114 | * coefIcoor = (2^30 -1)/D
|
---|
[5120] | 115 | */
|
---|
| 116 | coefIcoor=(MaxICoor)/(Max(pmax.x-pmin.x,pmax.y-pmin.y));
|
---|
[13036] | 117 | if(coefIcoor<=0) _error_("coefIcoor should be positive");
|
---|
[3913] | 118 | }
|
---|
| 119 | else{
|
---|
[13036] | 120 | _error_("No BamgVertex provided");
|
---|
[3913] | 121 | }
|
---|
| 122 |
|
---|
| 123 | //Edges
|
---|
| 124 | if (bamggeom->Edges){
|
---|
[5143] | 125 | R2 zerovector(0,0);
|
---|
| 126 | double* verticeslength=NULL;
|
---|
[3913] | 127 |
|
---|
[12509] | 128 | if(verbose>5) _printLine_(" processing Edges");
|
---|
[13036] | 129 | if (bamggeom->EdgesSize[1]!=3) _error_("Edges should have 3 columns");
|
---|
[5573] | 130 | edges = new GeomEdge[nbe];
|
---|
[3913] | 131 |
|
---|
[5143] | 132 | //initialize verticeslength (sum of the lengths of the edges holding vertex)
|
---|
| 133 | verticeslength = new double[nbv];
|
---|
| 134 | for(i=0;i<nbv;i++) verticeslength[i]=0;
|
---|
[3913] | 135 |
|
---|
[5143] | 136 | /*Loop over the edges*/
|
---|
[3913] | 137 | for (i=0;i<nbe;i++){
|
---|
| 138 |
|
---|
| 139 | i1=(int)bamggeom->Edges[i*3+0]-1; //-1 for C indexing
|
---|
| 140 | i2=(int)bamggeom->Edges[i*3+1]-1; //-1 for C indexing
|
---|
[5143] | 141 | edges[i].v[0]= vertices + i1; //pointer toward vertex i1 (=&vertices[i1])
|
---|
| 142 | edges[i].v[1]= vertices + i2; //pointer toward vertex i2
|
---|
[5148] | 143 | edges[i].ReferenceNumber=(long)bamggeom->Edges[i*3+2];
|
---|
[3913] | 144 |
|
---|
| 145 | //get length of edge
|
---|
[5143] | 146 | R2 x12=vertices[i2].r-vertices[i1].r;
|
---|
[3913] | 147 | double l12=sqrt((x12,x12));
|
---|
| 148 | Hmin=Min(Hmin,l12);
|
---|
| 149 |
|
---|
| 150 | //initialize other fields
|
---|
[5143] | 151 | edges[i].tg[0]=zerovector;
|
---|
| 152 | edges[i].tg[1]=zerovector;
|
---|
[5401] | 153 | edges[i].AdjVertexIndex[0] = edges[i].AdjVertexIndex[1] = -1;
|
---|
[5143] | 154 | edges[i].Adj[0] = edges[i].Adj[1] = NULL;
|
---|
[5146] | 155 | edges[i].type = 0;
|
---|
[3913] | 156 |
|
---|
| 157 | //Cracked?
|
---|
[5148] | 158 | if (edges[i].ReferenceNumber!=1) edges[i].SetCracked();
|
---|
[3913] | 159 |
|
---|
| 160 | //prepare metric
|
---|
| 161 | vertices[i1].color++;
|
---|
| 162 | vertices[i2].color++;
|
---|
[5143] | 163 | verticeslength[i1] += l12;
|
---|
| 164 | verticeslength[i2] += l12;
|
---|
[3913] | 165 | }
|
---|
| 166 |
|
---|
[5397] | 167 | // definition the default of the given mesh size
|
---|
[5143] | 168 | for (i=0;i<nbv;i++) {
|
---|
| 169 | if (vertices[i].color > 0)
|
---|
| 170 | vertices[i].m=Metric(verticeslength[i] /(double) vertices[i].color);
|
---|
| 171 | else
|
---|
| 172 | vertices[i].m=Metric(Hmin);
|
---|
| 173 | }
|
---|
| 174 | delete [] verticeslength;
|
---|
[3913] | 175 |
|
---|
| 176 | }
|
---|
| 177 | else{
|
---|
[13036] | 178 | _error_("No edges provided");
|
---|
[3913] | 179 | }
|
---|
| 180 |
|
---|
| 181 | //hVertices
|
---|
[5187] | 182 | if(bamgopts->hVertices && bamgopts->hVerticesSize[0]==nbv){
|
---|
[12509] | 183 | if(verbose>5) _printLine_(" processing hVertices");
|
---|
[3913] | 184 | for (i=0;i< nbv;i++){
|
---|
[13255] | 185 | if (!xIsNan<IssmPDouble>(bamgopts->hVertices[i])){
|
---|
[5187] | 186 | vertices[i].m=Metric((double)bamgopts->hVertices[i]);
|
---|
[3913] | 187 | }
|
---|
| 188 | }
|
---|
| 189 | }
|
---|
| 190 |
|
---|
| 191 | //MetricVertices
|
---|
[5187] | 192 | if(bamgopts->metric && bamgopts->metric[0]==nbv){
|
---|
[12509] | 193 | if(verbose>5) _printLine_(" processing MetricVertices");
|
---|
[3913] | 194 | for (i=0;i< nbv;i++) {
|
---|
[5187] | 195 | vertices[i].m = Metric((double)bamgopts->metric[i*3+0],(double)bamgopts->metric[i*3+1],(double)bamgopts->metric[i*3+2]);
|
---|
[3913] | 196 | }
|
---|
| 197 | }
|
---|
| 198 |
|
---|
[5091] | 199 | //MaxCornerAngle
|
---|
| 200 | if (bamgopts->MaxCornerAngle){
|
---|
[12509] | 201 | if(verbose>5) _printLine_(" processing MaxCornerAngle");
|
---|
[5091] | 202 | MaxCornerAngle=bamgopts->MaxCornerAngle*Pi/180;
|
---|
[3913] | 203 | }
|
---|
| 204 |
|
---|
| 205 | //TangentAtEdges
|
---|
| 206 | if (bamggeom->TangentAtEdges){
|
---|
[12509] | 207 | if(verbose>5) _printString_(" processing TangentAtEdges");
|
---|
[13036] | 208 | if (bamggeom->TangentAtEdgesSize[1]!=4) _error_("TangentAtEdges should have 4 columns");
|
---|
[3913] | 209 | int n,i,j,k;
|
---|
| 210 | R2 tg;
|
---|
| 211 |
|
---|
| 212 | n=bamggeom->TangentAtEdgesSize[0];
|
---|
| 213 | for (k=0;k<n;k++) {
|
---|
| 214 | i=(int)bamggeom->TangentAtEdges[k*4+0]-1; //for C indexing
|
---|
| 215 | j=(int)bamggeom->TangentAtEdges[k*4+1]-1; //for C indexing
|
---|
| 216 | tg.x=bamggeom->TangentAtEdges[k*4+2];
|
---|
| 217 | tg.y=bamggeom->TangentAtEdges[k*4+3];
|
---|
[13036] | 218 | if (i<0 || i>=nbe) _error_("TangentAtEdges first index exceeds matrix dimension");
|
---|
| 219 | if (j!=0 && j!=1) _error_("TangentAtEdges second index should be 1 or 2 only");
|
---|
[3913] | 220 | edges[i].tg[j] = tg;
|
---|
| 221 | }
|
---|
| 222 | }
|
---|
| 223 |
|
---|
| 224 | //Corners
|
---|
| 225 | if(bamggeom->Corners){
|
---|
[12509] | 226 | if(verbose>5) _printString_(" processing Corners");
|
---|
[13036] | 227 | if (bamggeom->CornersSize[1]!=1) _error_("Corners should have 1 column");
|
---|
[3913] | 228 | n=bamggeom->CornersSize[0];
|
---|
| 229 | for (i=0;i<n;i++) {
|
---|
| 230 | j=(int)bamggeom->Corners[i]-1; //for C indexing
|
---|
[13036] | 231 | if (j>nbv-1 || j<0) _error_("Bad corner definition: should in [0 " << nbv << "]");
|
---|
[5143] | 232 | /*Required => at the same time SetRequired and SetCorner*/
|
---|
[3913] | 233 | vertices[j].SetCorner();
|
---|
[5143] | 234 | vertices[j].SetRequired();
|
---|
| 235 | }
|
---|
[3913] | 236 | }
|
---|
| 237 |
|
---|
| 238 | //RequiredVertices
|
---|
| 239 | if(bamggeom->RequiredVertices){
|
---|
[12509] | 240 | if(verbose>5) _printLine_(" processing RequiredVertices");
|
---|
[13036] | 241 | if (bamggeom->RequiredVerticesSize[1]!=1) _error_("RequiredVertices should have 1 column");
|
---|
[3913] | 242 | n=bamggeom->RequiredVerticesSize[0];
|
---|
| 243 | for (i=0;i<n;i++) {
|
---|
| 244 | j=(int)bamggeom->RequiredVertices[i]-1; //for C indexing
|
---|
[13036] | 245 | if (j>nbv-1 || j<0) _error_("Bad RequiredVerticess definition: should in [0 " << nbv << "]");
|
---|
[5143] | 246 | vertices[j].SetRequired();
|
---|
| 247 | }
|
---|
[3913] | 248 | }
|
---|
| 249 |
|
---|
| 250 | //RequiredEdges
|
---|
| 251 | if(bamggeom->RequiredEdges){
|
---|
[12509] | 252 | if(verbose>5) _printLine_(" processing RequiredEdges");
|
---|
[13036] | 253 | if (bamggeom->RequiredEdgesSize[1]!=1) _error_("RequiredEdges should have 1 column");
|
---|
[3913] | 254 | n=bamggeom->RequiredEdgesSize[0];
|
---|
| 255 | for (i=0;i<n;i++) {
|
---|
| 256 | j=(int)bamggeom->RequiredEdges[i]-1; //for C indexing
|
---|
[13036] | 257 | if (j>nbe-1 || j<0) _error_("Bad RequiredEdges definition: should in [0 " << nbe << "]");
|
---|
[3913] | 258 | edges[j].SetRequired();
|
---|
| 259 | }
|
---|
| 260 | }
|
---|
| 261 |
|
---|
| 262 | //SubDomain
|
---|
| 263 | if(bamggeom->SubDomains){
|
---|
[12509] | 264 | if(verbose>5) _printLine_(" processing SubDomains");
|
---|
[13036] | 265 | if (bamggeom->SubDomainsSize[1]!=4) _error_("SubDomains should have 4 columns");
|
---|
[5148] | 266 | nbsubdomains=bamggeom->SubDomainsSize[0];
|
---|
[5573] | 267 | subdomains = new GeomSubDomain[nbsubdomains];
|
---|
[5148] | 268 | for (i=0;i<nbsubdomains;i++){
|
---|
[3913] | 269 | i0=(int)bamggeom->SubDomains[i*4+0];
|
---|
| 270 | i1=(int)bamggeom->SubDomains[i*4+1];
|
---|
| 271 | i2=(int)bamggeom->SubDomains[i*4+2];
|
---|
| 272 | i3=(int)bamggeom->SubDomains[i*4+3];
|
---|
[13036] | 273 | if (i0!=2) _error_("Bad Subdomain definition: first number should be 2 (for Edges)");
|
---|
| 274 | if (i1>nbe || i1<=0) _error_("Bad Subdomain definition: second number should in [1 " << nbe << "] (edge number)");
|
---|
[3913] | 275 | subdomains[i].edge=edges + (i1-1);
|
---|
[5148] | 276 | subdomains[i].direction = (int) i2;
|
---|
| 277 | subdomains[i].ReferenceNumber = i3;
|
---|
[3913] | 278 | }
|
---|
| 279 | }
|
---|
| 280 | }
|
---|
[12365] | 281 | /*}}}*/
|
---|
| 282 | /*FUNCTION Geometry::WriteGeometry{{{*/
|
---|
[3913] | 283 | void Geometry::WriteGeometry(BamgGeom* bamggeom, BamgOpts* bamgopts){
|
---|
| 284 |
|
---|
| 285 | int verbose;
|
---|
| 286 | int nbreq=0;
|
---|
| 287 | int nbreqv=0;
|
---|
| 288 | int nbtan=0;
|
---|
| 289 | int nbcracked=0;
|
---|
| 290 | int i,count;
|
---|
| 291 |
|
---|
| 292 | /*Get options*/
|
---|
| 293 | verbose=bamgopts->verbose;
|
---|
| 294 |
|
---|
| 295 | /*Vertices*/
|
---|
[12509] | 296 | if(verbose>5) _printLine_(" writing Vertices");
|
---|
[3913] | 297 | bamggeom->VerticesSize[0]=nbv;
|
---|
| 298 | bamggeom->VerticesSize[1]=3;
|
---|
| 299 | if (nbv){
|
---|
[12456] | 300 | bamggeom->Vertices=xNew<double>(3*nbv);
|
---|
[3913] | 301 | for (i=0;i<nbv;i++){
|
---|
| 302 | bamggeom->Vertices[i*3+0]=vertices[i].r.x;
|
---|
| 303 | bamggeom->Vertices[i*3+1]=vertices[i].r.y;
|
---|
[5148] | 304 | bamggeom->Vertices[i*3+2]=vertices[i].GetReferenceNumber();
|
---|
[3913] | 305 |
|
---|
| 306 | //update counters
|
---|
| 307 | if (vertices[i].Required()) nbreqv++;
|
---|
| 308 | }
|
---|
| 309 | }
|
---|
| 310 |
|
---|
| 311 | /*Edges*/
|
---|
[12509] | 312 | if(verbose>5) _printLine_(" writing Edges");
|
---|
[3913] | 313 | bamggeom->EdgesSize[0]=nbe;
|
---|
| 314 | bamggeom->EdgesSize[1]=3;
|
---|
| 315 | if (nbe){
|
---|
[12456] | 316 | bamggeom->Edges=xNew<double>(3*nbe);
|
---|
[3913] | 317 | for (i=0;i<nbe;i++){
|
---|
[5149] | 318 | bamggeom->Edges[i*3+0]=GetId(edges[i][0])+1; //back to Matlab indexing
|
---|
| 319 | bamggeom->Edges[i*3+1]=GetId(edges[i][1])+1; //back to Matlab indexing
|
---|
[5148] | 320 | bamggeom->Edges[i*3+2]=(double)edges[i].ReferenceNumber;
|
---|
[3913] | 321 |
|
---|
| 322 | //update counters
|
---|
| 323 | if (edges[i].Required()) nbreq++;
|
---|
| 324 | if (edges[i].TgA() && edges[i][0].Corner()) nbtan++;
|
---|
| 325 | if (edges[i].TgB() && edges[i][1].Corner()) nbtan++;
|
---|
| 326 | }
|
---|
| 327 | }
|
---|
| 328 |
|
---|
| 329 | /*RequiredEdges*/
|
---|
[12509] | 330 | if(verbose>5) _printLine_(" writing " << nbreq << " RequiredEdges");
|
---|
[3913] | 331 | bamggeom->RequiredEdgesSize[0]=nbreq;
|
---|
| 332 | bamggeom->RequiredEdgesSize[1]=1;
|
---|
| 333 | if (nbreq){
|
---|
[12456] | 334 | bamggeom->RequiredEdges=xNew<double>(1*nbreq);
|
---|
[3913] | 335 | count=0;
|
---|
| 336 | for (i=0;i<nbe;i++){
|
---|
| 337 | if (edges[i].Required()){
|
---|
| 338 | bamggeom->RequiredEdges[count]=i+1; //back to Matlab indexing
|
---|
| 339 | count=count+1;
|
---|
| 340 | }
|
---|
| 341 | }
|
---|
| 342 | }
|
---|
| 343 |
|
---|
| 344 | //No corners
|
---|
| 345 |
|
---|
| 346 | /*RequiredVertices*/
|
---|
[12509] | 347 | if(verbose>5) _printLine_(" writing " << nbreqv << " RequiredVertices");
|
---|
[3913] | 348 | bamggeom->RequiredVerticesSize[0]=nbreqv;
|
---|
| 349 | bamggeom->RequiredVerticesSize[1]=1;
|
---|
| 350 | if (nbreqv){
|
---|
[12456] | 351 | bamggeom->RequiredVertices=xNew<double>(1*nbreqv);
|
---|
[3913] | 352 | count=0;
|
---|
| 353 | for (i=0;i<nbv;i++){
|
---|
| 354 | if (vertices[i].Required()){
|
---|
| 355 | bamggeom->RequiredVertices[count]=i+1; //back to Matlab indexing
|
---|
| 356 | count=count+1;
|
---|
| 357 | }
|
---|
| 358 | }
|
---|
| 359 | }
|
---|
| 360 |
|
---|
| 361 | /*SubDomains*/
|
---|
[12509] | 362 | if(verbose>5) _printLine_(" writing SubDomains");
|
---|
[5148] | 363 | bamggeom->SubDomainsSize[0]=nbsubdomains;
|
---|
[3913] | 364 | bamggeom->SubDomainsSize[1]=4;
|
---|
[5148] | 365 | if (nbsubdomains){
|
---|
[12456] | 366 | bamggeom->SubDomains=xNew<double>(4*nbsubdomains);
|
---|
[5148] | 367 | for (i=0;i<nbsubdomains;i++){
|
---|
[3913] | 368 | bamggeom->SubDomains[4*i+0]=2;
|
---|
[5149] | 369 | bamggeom->SubDomains[4*i+1]=GetId(subdomains[i].edge)+1; //back to Matlab indexing
|
---|
[5148] | 370 | bamggeom->SubDomains[4*i+2]=subdomains[i].direction;
|
---|
| 371 | bamggeom->SubDomains[4*i+3]=subdomains[i].ReferenceNumber;
|
---|
[3913] | 372 | }
|
---|
| 373 | }
|
---|
| 374 |
|
---|
| 375 | /*TangentAtEdges*/
|
---|
[12509] | 376 | if(verbose>5) _printLine_(" writing TangentAtEdges");
|
---|
[3913] | 377 | bamggeom->TangentAtEdgesSize[0]=nbtan;
|
---|
| 378 | bamggeom->TangentAtEdgesSize[1]=4;
|
---|
| 379 | if (nbtan){
|
---|
[12456] | 380 | bamggeom->TangentAtEdges=xNew<double>(4*nbtan);
|
---|
[3913] | 381 | count=0;
|
---|
| 382 | for (i=0;i<nbe;i++){
|
---|
| 383 | if (edges[i].TgA() && edges[i][0].Corner()){
|
---|
| 384 | bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
|
---|
| 385 | bamggeom->TangentAtEdges[4*i+1]=1;
|
---|
| 386 | bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[0].x;
|
---|
| 387 | bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[0].y;
|
---|
| 388 | }
|
---|
| 389 | if (edges[i].TgB() && edges[i][1].Corner()){
|
---|
| 390 | bamggeom->TangentAtEdges[4*i+0]=i+1; //back to Matlab indexing
|
---|
| 391 | bamggeom->TangentAtEdges[4*i+1]=2;
|
---|
| 392 | bamggeom->TangentAtEdges[4*i+2]=edges[i].tg[1].x;
|
---|
[5016] | 393 | bamggeom->TangentAtEdges[4*i+3]=edges[i].tg[1].y;
|
---|
[3913] | 394 | }
|
---|
| 395 | count=count+1;
|
---|
| 396 | }
|
---|
| 397 | }
|
---|
| 398 | }
|
---|
[12365] | 399 | /*}}}*/
|
---|
[3913] | 400 |
|
---|
| 401 | /*Methods*/
|
---|
[12365] | 402 | /*FUNCTION Geometry::Echo {{{*/
|
---|
[5400] | 403 | void Geometry::Echo(void){
|
---|
| 404 |
|
---|
[12509] | 405 | _printLine_("Geometry:");
|
---|
| 406 | _printLine_(" nbv (number of vertices) : " << nbv);
|
---|
| 407 | _printLine_(" nbe (number of edges) : " << nbe);
|
---|
| 408 | _printLine_(" nbsubdomains: " << nbsubdomains);
|
---|
| 409 | _printLine_(" nbcurves: " << nbcurves);
|
---|
| 410 | _printLine_(" vertices: " << vertices);
|
---|
| 411 | _printLine_(" edges: " << edges);
|
---|
| 412 | _printLine_(" quadtree: " << quadtree);
|
---|
| 413 | _printLine_(" subdomains: " << subdomains);
|
---|
| 414 | _printLine_(" curves: " << curves);
|
---|
| 415 | _printLine_(" pmin (x,y): (" << pmin.x << " " << pmin.y << ")");
|
---|
| 416 | _printLine_(" pmax (x,y): (" << pmax.x << " " << pmax.y << ")");
|
---|
| 417 | _printLine_(" coefIcoor: " << coefIcoor);
|
---|
| 418 | _printLine_(" MaxCornerAngle: " << MaxCornerAngle);
|
---|
[5400] | 419 |
|
---|
| 420 | return;
|
---|
| 421 | }
|
---|
| 422 | /*}}}*/
|
---|
[12365] | 423 | /*FUNCTION Geometry::Init{{{*/
|
---|
[5400] | 424 | void Geometry::Init(void){
|
---|
| 425 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/EmptyGeometry)*/
|
---|
| 426 |
|
---|
| 427 | NbRef=0;
|
---|
| 428 | nbv=0;
|
---|
| 429 | nbe=0;
|
---|
| 430 | quadtree=NULL;
|
---|
| 431 | curves=NULL;
|
---|
| 432 | edges=NULL;
|
---|
| 433 | vertices=NULL;
|
---|
| 434 | nbsubdomains=0;
|
---|
| 435 | nbcurves=0;
|
---|
| 436 | subdomains=NULL;
|
---|
| 437 | MaxCornerAngle = 10*Pi/180; //default is 10 degres
|
---|
| 438 | }
|
---|
[12365] | 439 | /*}}}*/
|
---|
| 440 | /*FUNCTION Geometry::MinimalHmin{{{*/
|
---|
[5400] | 441 | double Geometry::MinimalHmin() {
|
---|
| 442 | /* coeffIcoor = (2^30-1)/D
|
---|
| 443 | * We cannot go beyond hmin = D/2^30 because of the quadtree
|
---|
| 444 | * hmin is therefore approximately 2/coeffIcoor */
|
---|
| 445 | return 2.0/coefIcoor;
|
---|
| 446 | }/*}}}*/
|
---|
[12365] | 447 | /*FUNCTION Geometry::MaximalHmax{{{*/
|
---|
[5400] | 448 | double Geometry::MaximalHmax() {
|
---|
| 449 | return Max(pmax.x-pmin.x,pmax.y-pmin.y);
|
---|
| 450 | }/*}}}*/
|
---|
[12365] | 451 | /*FUNCTION Geometry::GetId(const GeomVertex &t){{{*/
|
---|
[5573] | 452 | long Geometry::GetId(const GeomVertex & t) const {
|
---|
[5400] | 453 | return &t - vertices;
|
---|
| 454 | }/*}}}*/
|
---|
[12365] | 455 | /*FUNCTION Geometry::GetId(const GeomVertex * t){{{*/
|
---|
[5573] | 456 | long Geometry::GetId(const GeomVertex * t) const {
|
---|
[5400] | 457 | return t - vertices;
|
---|
| 458 | }/*}}}*/
|
---|
[12365] | 459 | /*FUNCTION Geometry::GetId(const GeomEdge & t){{{*/
|
---|
[5573] | 460 | long Geometry::GetId(const GeomEdge & t) const {
|
---|
[5400] | 461 | return &t - edges;
|
---|
| 462 | }/*}}}*/
|
---|
[12365] | 463 | /*FUNCTION Geometry::GetId(const GeomEdge * t){{{*/
|
---|
[5573] | 464 | long Geometry::GetId(const GeomEdge * t) const {
|
---|
[5400] | 465 | return t - edges;
|
---|
| 466 | }/*}}}*/
|
---|
[12365] | 467 | /*FUNCTION Geometry::GetId(const Curve * c){{{*/
|
---|
[5400] | 468 | long Geometry::GetId(const Curve * c) const {
|
---|
| 469 | return c - curves;
|
---|
| 470 | }/*}}}*/
|
---|
[12365] | 471 | /*FUNCTION Geometry::Containing{{{*/
|
---|
[5573] | 472 | GeomEdge* Geometry::Containing(const R2 P, GeomEdge * start) const {
|
---|
[5400] | 473 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/Contening)*/
|
---|
| 474 |
|
---|
[5573] | 475 | GeomEdge* on =start,* pon=0;
|
---|
[5400] | 476 | // walk with the cos on geometry
|
---|
| 477 | int counter=0;
|
---|
| 478 | while(pon != on){
|
---|
| 479 | counter++;
|
---|
[6412] | 480 | _assert_(counter<100);
|
---|
[5400] | 481 | pon = on;
|
---|
| 482 | R2 A= (*on)[0];
|
---|
| 483 | R2 B= (*on)[1];
|
---|
| 484 | R2 AB = B-A;
|
---|
| 485 | R2 AP = P-A;
|
---|
| 486 | R2 BP = P-B;
|
---|
| 487 | if ( (AB,AP) < 0)
|
---|
| 488 | on = on->Adj[0];
|
---|
| 489 | else if ( (AB,BP) > 0)
|
---|
| 490 | on = on->Adj[1];
|
---|
| 491 | else
|
---|
| 492 | return on;
|
---|
| 493 | }
|
---|
| 494 | return on;
|
---|
| 495 | }
|
---|
[12365] | 496 | /*}}}*/
|
---|
| 497 | /*FUNCTION Geometry::PostRead{{{*/
|
---|
[5397] | 498 | void Geometry::PostRead(){
|
---|
[3913] | 499 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/AfterRead)*/
|
---|
| 500 |
|
---|
[5177] | 501 | long i,j,k;
|
---|
| 502 | int jj;
|
---|
| 503 | long *head_v = new long[nbv];
|
---|
| 504 | long *next_p = new long[2*nbe];
|
---|
| 505 | float *eangle = new float[nbe];
|
---|
| 506 | double eps = 1e-20;
|
---|
[12218] | 507 | BamgQuadtree quadtree; // build quadtree to find duplicates
|
---|
[5177] | 508 | BamgVertex *v0 = vertices;
|
---|
[5573] | 509 | GeomVertex *v0g = (GeomVertex*) (void*)v0;
|
---|
[3913] | 510 |
|
---|
| 511 | k=0;
|
---|
| 512 |
|
---|
| 513 | //build quadtree for this geometry
|
---|
| 514 | for (i=0;i<nbv;i++){
|
---|
| 515 |
|
---|
| 516 | /*build integer coordinates (non unique)
|
---|
[5400] | 517 | these coordinates are used by the quadtree to group
|
---|
| 518 | the vertices by groups of 5:
|
---|
| 519 | All the coordinates are transformed to ]0,1[^2
|
---|
| 520 | then, the integer coordinates are computed using
|
---|
| 521 | the transformation ]0,1[^2 -> [0 2^30-1[^2 for a quadtree of depth 30*/
|
---|
[5180] | 522 | vertices[i].i=R2ToI2(vertices[i].r);
|
---|
[3913] | 523 |
|
---|
[5177] | 524 | /*find nearest vertex already present in the quadtree (NULL if empty)*/
|
---|
[5120] | 525 | BamgVertex* v=quadtree.NearestVertex(vertices[i].i.x,vertices[i].i.y);
|
---|
[3913] | 526 |
|
---|
[5177] | 527 | /*if there is a vertex found that is to close to vertices[i] -> error*/
|
---|
| 528 | if( v && Norme1(v->r - vertices[i].r) < eps ){
|
---|
[12509] | 529 | _printLine_("reference numbers: " << v->ReferenceNumber << " " << vertices[i].ReferenceNumber);
|
---|
| 530 | _printLine_("Id: " << i+1);
|
---|
[13385] | 531 | _printLine_("Coords: ["<<v->r.x<<" "<<v->r.y<<"] ["<<vertices[i].r.x<<" "<<vertices[i].r.y<<"]");
|
---|
| 532 |
|
---|
[9306] | 533 | delete [] next_p;
|
---|
| 534 | delete [] head_v;
|
---|
| 535 | delete [] eangle;
|
---|
[13036] | 536 | _error_("two points of the geometry are very closed to each other (see reference numbers above)");
|
---|
[3913] | 537 | }
|
---|
| 538 |
|
---|
[5177] | 539 | /*Add vertices[i] to the quadtree*/
|
---|
[5124] | 540 | quadtree.Add(vertices[i]);
|
---|
[3913] | 541 | }
|
---|
| 542 |
|
---|
| 543 | /* Here we use a powerful chaining algorithm
|
---|
| 544 | *
|
---|
| 545 | * 1. What is a chaining algorithm?
|
---|
| 546 | *
|
---|
| 547 | * If F is a function that goes from i in [0 n] to j in [0 m]
|
---|
| 548 | * and we want to compute the reciprocal function F-1 of F
|
---|
| 549 | * (what are the antecedents of a given j in Im(F) )
|
---|
| 550 | * We use 2 lists:
|
---|
| 551 | * head_F[j] that holds the head of lists
|
---|
| 552 | * next_F[i] that holds the list of elements that have the same image
|
---|
| 553 | *
|
---|
| 554 | * Example:
|
---|
| 555 | * i1, i2, ..., ip in [0,n] are all antecedents of a given j in [0 m]
|
---|
| 556 | * head_F[j] = ip
|
---|
| 557 | * next_F[ip]= ip-1
|
---|
| 558 | * ....
|
---|
| 559 | * next_F[i2]= i1
|
---|
| 560 | * next_F[i1]= -1 //end of the list
|
---|
| 561 | *
|
---|
| 562 | * Algorithm:
|
---|
| 563 | * for(j=0;j<m;j++) head_F[j] = -1 //initialization
|
---|
| 564 | * for(i=0;i<n;i++){
|
---|
| 565 | * j=F[i];
|
---|
| 566 | * next_F[i]= head_F[j];
|
---|
| 567 | * head_F[j]=i;
|
---|
| 568 | * }
|
---|
| 569 | *
|
---|
| 570 | * Then, we can go through all the elements that have for image j:
|
---|
| 571 | * for(i=head_F[j]; i!=-1; i=next_F[i])
|
---|
| 572 | * initialization of i by i=head_F[j]
|
---|
| 573 | * stop the loop when i=-1 (end of the chain)
|
---|
| 574 | * iterate using i=next_F[i] (next element that have for image j)
|
---|
| 575 | *
|
---|
| 576 | * 2. How to use this algorithm here?
|
---|
| 577 | *
|
---|
| 578 | * Here F is a function that associates two vertices v0 and v1 for a given edge E
|
---|
| 579 | * We want to build the reciprocal function: what are the edges that contains
|
---|
| 580 | * a vertex v?
|
---|
| 581 | * To do so, we use the same chaining algorithm but there is a difficulty
|
---|
[5339] | 582 | * coming from the fact that for F we have a couple of vertices and not one
|
---|
| 583 | * vertex.
|
---|
| 584 | * To overcome this difficulty, we use a global indexing exactly like in
|
---|
[3913] | 585 | * C/C++ so that
|
---|
| 586 | * a member of a 2-column-table can be described by one index p=i*2+j
|
---|
| 587 | * i=(int)p/2 line number of p
|
---|
[5339] | 588 | * j=p%2 column number of p
|
---|
[3913] | 589 | *
|
---|
| 590 | * Algorithm:
|
---|
| 591 | * for(i=0;i<nbv;i++) head_v[i] = -1 //initialization
|
---|
| 592 | * for(i=0;i<nbe;i++){
|
---|
| 593 | * for(j=0;j<2;j++){
|
---|
| 594 | * p=2*i+j;
|
---|
| 595 | * v=edges(i,j);
|
---|
| 596 | * next_p[p]= head_v[v];
|
---|
| 597 | * head_v[v]=p;
|
---|
| 598 | * }
|
---|
| 599 | * }
|
---|
| 600 | */
|
---|
| 601 |
|
---|
| 602 | //initialize head_v as -1
|
---|
| 603 | for (i=0;i<nbv;i++) head_v[i]=-1;
|
---|
| 604 | k=0;
|
---|
| 605 | for (i=0;i<nbe;i++) {
|
---|
| 606 | //compute vector of edge i that goes from vertex 0 to vertex 1
|
---|
| 607 | R2 v10=edges[i].v[1]->r - edges[i].v[0]->r;
|
---|
| 608 | double lv10=Norme2(v10);
|
---|
| 609 | //check that its length is not 0
|
---|
[9309] | 610 | if(lv10==0){
|
---|
| 611 | delete [] next_p;
|
---|
| 612 | delete [] head_v;
|
---|
| 613 | delete [] eangle;
|
---|
[13036] | 614 | _error_("Length of edge " << i << " is 0");
|
---|
[9309] | 615 | }
|
---|
[3913] | 616 | //compute angle in [-Pi Pi]
|
---|
| 617 | eangle[i] = atan2(v10.y,v10.x);
|
---|
| 618 | //build chains head_v and next_p
|
---|
| 619 | for (j=0;j<2;j++){
|
---|
[5149] | 620 | long v=GetId(edges[i].v[j]);
|
---|
[3913] | 621 | next_p[k]=head_v[v];
|
---|
| 622 | head_v[v]=k++; //post increment: head_v[v]=k; and then k=k+1;
|
---|
| 623 | }
|
---|
| 624 | }
|
---|
| 625 |
|
---|
| 626 | //sort head_v by order of increasing edges angle
|
---|
| 627 | for (i=0;i<nbv;i++) {
|
---|
[5339] | 628 | int exch=1,ord=0;
|
---|
[3913] | 629 |
|
---|
| 630 | //exchange vertices position in head_v and next_p till tey are sorted
|
---|
| 631 | while (exch){
|
---|
[5339] | 632 | long *p=head_v+i;
|
---|
| 633 | long *po=p;
|
---|
| 634 | long n=*p;
|
---|
[3913] | 635 | register float angleold=-1000 ; // angle = - infinity
|
---|
| 636 | ord=0; exch=0;
|
---|
| 637 |
|
---|
[5339] | 638 | // loop over the edges that contain the vertex i (till -1)
|
---|
[3913] | 639 | while (n >=0){
|
---|
| 640 | ord++;
|
---|
[5339] | 641 | long i1=n/2; // i1 = floor (n/2) -> Edge number
|
---|
| 642 | long j1=n%2; // j1 = 1 if n is odd -> Vertex index for this edge (0 or 1)
|
---|
| 643 | long* pn=next_p+n;
|
---|
[3913] | 644 |
|
---|
[5339] | 645 | //Next vertex index
|
---|
[3913] | 646 | n=*pn;
|
---|
| 647 |
|
---|
| 648 | //compute angle between horizontal axis and v0->v1
|
---|
| 649 | float angle = j1 ? OppositeAngle(eangle[i1]): eangle[i1];
|
---|
| 650 |
|
---|
| 651 | //exchange if the current edge angle is smaller than the previous one
|
---|
| 652 | if (angleold > angle){
|
---|
| 653 | exch=1;
|
---|
| 654 | *pn=*po; // next_p[n] = n + 1
|
---|
| 655 | *po=*p; //
|
---|
| 656 | *p=n; // next_p[n+1] = n
|
---|
| 657 | po=pn; // po now point toward pn (invert next and current)
|
---|
| 658 | }
|
---|
| 659 |
|
---|
| 660 | //else, continue going to the next edge position
|
---|
| 661 | else{ // to have : po -> p -> pn
|
---|
| 662 | angleold=angle; // update maximum angle
|
---|
| 663 | po=p; // po now point toward p (current position)
|
---|
| 664 | p=pn; // p now point toward pn (next position)
|
---|
| 665 | }
|
---|
| 666 | }
|
---|
| 667 | }
|
---|
| 668 |
|
---|
[5339] | 669 | // angular test on current vertex to guess whether it is a corner (ord = number of edges holding i)
|
---|
| 670 | if(ord==2) {
|
---|
[3913] | 671 | long n1 = head_v[i];
|
---|
| 672 | long n2 = next_p[n1];
|
---|
| 673 | long i1 = n1/2, i2 = n2/2; // edge number
|
---|
| 674 | long j1 = n1%2, j2 = n2%2; // vertex in the edge
|
---|
| 675 | float angle1= j1 ? OppositeAngle(eangle[i1]) : eangle[i1];
|
---|
| 676 | float angle2= !j2 ? OppositeAngle(eangle[i2]) : eangle[i2];
|
---|
| 677 | float da12 = Abs(angle2-angle1);
|
---|
[5091] | 678 | if (( da12 >= MaxCornerAngle ) && (da12 <= 2*Pi -MaxCornerAngle)) {
|
---|
[3913] | 679 | vertices[i].SetCorner() ;
|
---|
| 680 | }
|
---|
[5339] | 681 | // if the edge type/referencenumber a changing then is SetRequired();
|
---|
[5146] | 682 | if (edges[i1].type != edges[i2].type || edges[i1].Required()){
|
---|
[3913] | 683 | vertices[i].SetRequired();
|
---|
| 684 | }
|
---|
[5148] | 685 | if (edges[i1].ReferenceNumber != edges[i2].ReferenceNumber) {
|
---|
[3913] | 686 | vertices[i].SetRequired();
|
---|
| 687 | }
|
---|
| 688 | }
|
---|
| 689 | if(ord != 2) {
|
---|
| 690 | vertices[i].SetCorner();
|
---|
| 691 | }
|
---|
| 692 |
|
---|
[5339] | 693 | /*close the list around the vertex to have a circular loop*/
|
---|
[3913] | 694 | long no=-1, ne = head_v[i];
|
---|
| 695 | while (ne >=0) ne = next_p[no=ne];
|
---|
| 696 | if(no>=0) next_p[no] = head_v[i];
|
---|
| 697 | }
|
---|
| 698 |
|
---|
[5339] | 699 | /*Check that the list makes sense (we have all the time the same vertex)
|
---|
| 700 | * and build adjacent edges*/
|
---|
[3913] | 701 | k =0;
|
---|
| 702 | for (i=0;i<nbe;i++){
|
---|
| 703 | for (j=0;j<2;j++){
|
---|
[5339] | 704 |
|
---|
[3913] | 705 | long n1 = next_p[k++];
|
---|
| 706 | long i1 = n1/2 ,j1=n1%2;
|
---|
[5339] | 707 |
|
---|
[13036] | 708 | if( edges[i1].v[j1] != edges[i].v[j]) _error_("Problem while processing edges: check the edge list");
|
---|
[5339] | 709 |
|
---|
[3913] | 710 | edges[i1].Adj[j1] = edges + i;
|
---|
[5401] | 711 | edges[i1].AdjVertexIndex[j1] = j;
|
---|
[3913] | 712 | }
|
---|
| 713 | }
|
---|
| 714 |
|
---|
[5397] | 715 | /* generation of all the tangents*/
|
---|
[3913] | 716 | for (i=0;i<nbe;i++) {
|
---|
[5339] | 717 | R2 AB =edges[i].v[1]->r -edges[i].v[0]->r;// AB = vertex0 -> vertex1
|
---|
| 718 | double lAB=Norme2(AB); // Get length of AB
|
---|
| 719 | double ltg2[2]={0.0}; // initialize tangent
|
---|
[3913] | 720 |
|
---|
| 721 | //loop over the 2 vertices of the edge
|
---|
[5339] | 722 | for (j=0;j<2;j++) {
|
---|
| 723 | R2 tg =edges[i].tg[j];
|
---|
[3913] | 724 | double ltg=Norme2(tg);
|
---|
| 725 |
|
---|
| 726 | //by default, tangent=[0 0]
|
---|
[5339] | 727 | if(ltg==0){
|
---|
[3913] | 728 | //if the current vertex of the edge is not a corner
|
---|
[5339] | 729 | if(!edges[i].v[j]->Corner()){
|
---|
| 730 | /*The tangent is set as the vector between the
|
---|
| 731 | * previous and next vertices connected to current vertex
|
---|
| 732 | * normed by the edge length*/
|
---|
[5401] | 733 | tg = edges[i].v[1-j]->r - edges[i].Adj[j]->v[1-edges[i].AdjVertexIndex[j]]->r;
|
---|
[3913] | 734 | ltg= Norme2(tg);
|
---|
| 735 | tg = tg *(lAB/ltg);
|
---|
| 736 | ltg= lAB;
|
---|
| 737 | }
|
---|
[5339] | 738 | //else: a Corner no tangent => nothing to do
|
---|
[3913] | 739 | }
|
---|
| 740 | else{
|
---|
| 741 | //tangent has already been computed
|
---|
| 742 | tg = tg*(lAB/ltg),ltg=lAB;
|
---|
| 743 | }
|
---|
| 744 |
|
---|
[5339] | 745 | ltg2[j] = ltg;
|
---|
[3913] | 746 |
|
---|
[5339] | 747 | if ((tg,AB)<0) tg = -tg;
|
---|
[3913] | 748 |
|
---|
[5339] | 749 | edges[i].tg[j]=tg;
|
---|
[3913] | 750 | }
|
---|
| 751 | if (ltg2[0]!=0) edges[i].SetTgA();
|
---|
| 752 | if (ltg2[1]!=0) edges[i].SetTgB();
|
---|
| 753 | }
|
---|
| 754 |
|
---|
[5397] | 755 | /* generation of all curves (from corner to corner)*/
|
---|
| 756 | /*We proceed in 2 steps: first allocate, second build*/
|
---|
[3913] | 757 | for (int step=0;step<2;step++){
|
---|
[5339] | 758 |
|
---|
[5397] | 759 | //unmark all edges
|
---|
[3913] | 760 | for (i=0;i<nbe;i++) edges[i].SetUnMark();
|
---|
[5397] | 761 | long nb_marked_edges=0;
|
---|
[5339] | 762 |
|
---|
[5397] | 763 | //initialize number of curves
|
---|
[5148] | 764 | nbcurves = 0;
|
---|
[4997] | 765 |
|
---|
[5397] | 766 | for (int level=0;level<2 && nb_marked_edges!=nbe;level++){
|
---|
[5339] | 767 | for (i=0;i<nbe;i++){
|
---|
[5397] | 768 |
|
---|
[5573] | 769 | GeomEdge & ei=edges[i];
|
---|
[5339] | 770 | for(j=0;j<2;j++){
|
---|
[5397] | 771 | /*If current edge ei is unmarked and (level=1 or vertex i is required (corner)):
|
---|
| 772 | * we do have the first edge of a new curve*/
|
---|
[5339] | 773 | if (!ei.Mark() && (level || ei[j].Required())) {
|
---|
| 774 | int k0=j,k1;
|
---|
[5573] | 775 | GeomEdge *e=&ei;
|
---|
| 776 | GeomVertex *a=(*e)(k0); // begin
|
---|
[5339] | 777 | if(curves){
|
---|
[5400] | 778 | curves[nbcurves].FirstEdge=e;
|
---|
| 779 | curves[nbcurves].FirstVertexIndex=k0;
|
---|
[5339] | 780 | }
|
---|
| 781 | int nee=0;
|
---|
[5397] | 782 | for(;;){
|
---|
[5339] | 783 | nee++;
|
---|
| 784 | k1 = 1-k0; // next vertex of the edge
|
---|
| 785 | e->SetMark();
|
---|
[5397] | 786 | nb_marked_edges++;
|
---|
[5339] | 787 | e->CurveNumber=nbcurves;
|
---|
[5573] | 788 | GeomVertex *b=(*e)(k1);
|
---|
[5397] | 789 |
|
---|
| 790 | //break if we have reached the other end of the curve
|
---|
| 791 | if (a==b || b->Required()){
|
---|
[5400] | 792 | if(curves){
|
---|
| 793 | curves[nbcurves].LastEdge=e;
|
---|
| 794 | curves[nbcurves].LastVertexIndex=k1;
|
---|
| 795 | }
|
---|
[5397] | 796 | break;
|
---|
| 797 | }
|
---|
| 798 | //else: go to next edge (adjacent)
|
---|
| 799 | else{
|
---|
[5401] | 800 | k0 = e->AdjVertexIndex[k1];// vertex in next edge
|
---|
[5397] | 801 | e = e->Adj[k1]; // next edge
|
---|
| 802 | }
|
---|
[5339] | 803 | }
|
---|
| 804 | nbcurves++;
|
---|
| 805 | if(level) a->SetRequired();
|
---|
| 806 | }
|
---|
| 807 | }
|
---|
| 808 | }
|
---|
| 809 | }
|
---|
[6412] | 810 | _assert_(nb_marked_edges && nbe);
|
---|
[5397] | 811 | //allocate if first step
|
---|
| 812 | if(step==0) curves=new Curve[nbcurves];
|
---|
[3913] | 813 | }
|
---|
| 814 |
|
---|
| 815 | /*clean up*/
|
---|
[4997] | 816 | delete [] next_p;
|
---|
| 817 | delete [] head_v;
|
---|
| 818 | delete [] eangle;
|
---|
[3913] | 819 |
|
---|
| 820 | }
|
---|
[12365] | 821 | /*}}}*/
|
---|
| 822 | /*FUNCTION Geometry::ProjectOnCurve {{{*/
|
---|
[5573] | 823 | GeomEdge* Geometry::ProjectOnCurve(const Edge &e,double s,BamgVertex &V,VertexOnGeom &GV) const {
|
---|
[3913] | 824 | /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, MeshGeom.cpp/ProjectOnCurve)*/
|
---|
| 825 | /*Add a vertex on an existing geometrical edge according to the metrics of the two vertices constituting the edge*/
|
---|
| 826 |
|
---|
| 827 | double save_s=s;
|
---|
| 828 | int NbTry=0;
|
---|
| 829 |
|
---|
| 830 | retry:
|
---|
| 831 |
|
---|
| 832 | s=save_s;
|
---|
[5573] | 833 | GeomEdge* on=e.GeomEdgeHook;
|
---|
[3913] | 834 | if (!on){
|
---|
[13036] | 835 | _error_("ProjectOnCurve error message: edge provided should be on geometry");
|
---|
[3913] | 836 | }
|
---|
[5573] | 837 | if (!e[0].GeomEdgeHook || !e[1].GeomEdgeHook){
|
---|
[13036] | 838 | _error_("ProjectOnCurve error message: at least one of the vertex of the edge provided is not on geometry");
|
---|
[3913] | 839 | }
|
---|
| 840 |
|
---|
| 841 | //Get the two vertices of the edge
|
---|
[5120] | 842 | const BamgVertex &v0=e[0];
|
---|
| 843 | const BamgVertex &v1=e[1];
|
---|
[3913] | 844 |
|
---|
| 845 | //Get position of V0, V1 and vector v0->v1
|
---|
| 846 | R2 V0=v0,V1=v1,V01=V1-V0;
|
---|
| 847 |
|
---|
| 848 | //Get geometrical vertices corresponding to v0 and v1
|
---|
[5573] | 849 | VertexOnGeom vg0=*v0.GeomEdgeHook, vg1=*v1.GeomEdgeHook;
|
---|
[3913] | 850 |
|
---|
| 851 | //build two pointers towrad current geometrical edge
|
---|
[5573] | 852 | GeomEdge *eg0=on, *eg1=on;
|
---|
[3913] | 853 |
|
---|
| 854 | //Get edge direction and swap v0 and v1 if necessary
|
---|
| 855 | R2 Ag=(R2)(*on)[0],Bg=(R2)(*on)[1],AB=Bg-Ag;
|
---|
| 856 | int OppositeSens = (V01,AB)<0;
|
---|
[5148] | 857 | int direction0=0,direction1=1;
|
---|
[3913] | 858 | if (OppositeSens) s=1-s,Exchange(vg0,vg1),Exchange(V0,V1);
|
---|
| 859 |
|
---|
| 860 | //Compute metric of new vertex as a linear interpolation of the two others
|
---|
| 861 | V.m=Metric(1.0-s,v0,s,v1);
|
---|
| 862 |
|
---|
| 863 | const int mxe=100;
|
---|
[5573] | 864 | GeomEdge* ge[mxe+1];
|
---|
[5148] | 865 | int directionge[mxe+1];
|
---|
[3913] | 866 | double lge[mxe+1];
|
---|
| 867 | int bge=mxe/2,tge=bge;
|
---|
[5573] | 868 | ge[bge] = e.GeomEdgeHook;
|
---|
[5148] | 869 | directionge[bge]=1;
|
---|
[3913] | 870 |
|
---|
[5573] | 871 | while (eg0!=(GeomEdge*)vg0 && (*eg0)(direction0)!=(GeomVertex*)vg0){
|
---|
[3913] | 872 | if (bge<=0) {
|
---|
| 873 | if(NbTry) {
|
---|
[12509] | 874 | _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
|
---|
| 875 | _printLine_("That bug might come from:");
|
---|
| 876 | _printLine_(" 1) a mesh edge containing more than " << mxe/2 << " geometrical edges");
|
---|
| 877 | _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
|
---|
| 878 | _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
|
---|
[13036] | 879 | _error_("see above");
|
---|
[5130] | 880 | }
|
---|
[3913] | 881 | NbTry++;
|
---|
| 882 | goto retry;
|
---|
| 883 | }
|
---|
[5573] | 884 | GeomEdge* tmpge = eg0;
|
---|
[5148] | 885 | ge[--bge] =eg0 = eg0->Adj[direction0];
|
---|
[6412] | 886 | _assert_(bge>=0 && bge<=mxe);
|
---|
[5401] | 887 | direction0 = 1-( directionge[bge] = tmpge->AdjVertexIndex[direction0]);
|
---|
[3913] | 888 | }
|
---|
[5573] | 889 | while (eg1 != (GeomEdge*) vg1 && (*eg1)(direction1) != (GeomVertex*) vg1) {
|
---|
[3913] | 890 | if(tge>=mxe ) {
|
---|
[12509] | 891 | _printLine_("WARNING: on the class Mesh before call Geometry::ProjectOnCurve is having issues (isn't it Eric?)");
|
---|
[3913] | 892 | NbTry++;
|
---|
| 893 | if (NbTry<2) goto retry;
|
---|
[12509] | 894 | _printLine_("Fatal Error: on the class Mesh before call Geometry::ProjectOnCurve");
|
---|
| 895 | _printLine_("That bug might come from:");
|
---|
| 896 | _printLine_(" 1) a mesh edge contening more than " << mxe/2 << " geometrical edges");
|
---|
| 897 | _printLine_(" 2) code bug : be sure that we call Mesh::SetVertexFieldOn() before");
|
---|
| 898 | _printLine_("To solve the problem do a coarsening of the geometrical mesh or change the constant value of mxe (dangerous)");
|
---|
[13036] | 899 | _error_("see above");
|
---|
[3913] | 900 | }
|
---|
[5573] | 901 | GeomEdge* tmpge = eg1;
|
---|
[5148] | 902 | ge[++tge] =eg1 = eg1->Adj[direction1];
|
---|
[5401] | 903 | directionge[tge]= direction1 = 1-tmpge->AdjVertexIndex[direction1];
|
---|
[6412] | 904 | _assert_(tge>=0 && tge<=mxe);
|
---|
[3913] | 905 | }
|
---|
| 906 |
|
---|
| 907 |
|
---|
[5573] | 908 | if ((*eg0)(direction0)==(GeomVertex*)vg0)
|
---|
[5148] | 909 | vg0=VertexOnGeom(*(BamgVertex*) vg0,*eg0,direction0); //vg0 = absisce
|
---|
[3913] | 910 |
|
---|
[5573] | 911 | if ((*eg1)(direction1)==(GeomVertex*)vg1)
|
---|
[5148] | 912 | vg1=VertexOnGeom(*(BamgVertex*) vg1,*eg1,direction1);
|
---|
[3913] | 913 |
|
---|
| 914 | double sg;
|
---|
| 915 | if (eg0 == eg1) {
|
---|
| 916 | register double s0=vg0,s1=vg1;
|
---|
| 917 | sg = s0*(1.0-s) + s*s1;
|
---|
| 918 | on=eg0;
|
---|
| 919 | }
|
---|
| 920 | else {
|
---|
| 921 | R2 AA=V0,BB;
|
---|
| 922 | double s0,s1;
|
---|
| 923 | int i=bge;
|
---|
| 924 | double ll=0;
|
---|
| 925 | for(i=bge;i<tge;i++){
|
---|
[6412] | 926 | _assert_(i>=0 && i<=mxe);
|
---|
[5148] | 927 | BB = (*ge[i])[directionge[i]];
|
---|
[3913] | 928 | lge[i]=ll += Norme2(AA-BB);
|
---|
| 929 | AA=BB ;}
|
---|
| 930 | lge[tge]=ll+=Norme2(AA-V1);
|
---|
| 931 | // search the geometrical edge
|
---|
[6412] | 932 | _assert_(s<=1.0);
|
---|
[3913] | 933 | double ls= s*ll;
|
---|
| 934 | on =0;
|
---|
| 935 | s0 = vg0;
|
---|
[5148] | 936 | s1= directionge[bge];
|
---|
[3913] | 937 | double l0=0,l1;
|
---|
| 938 | i=bge;
|
---|
| 939 | while ( (l1=lge[i]) < ls ) {
|
---|
[6412] | 940 | _assert_(i>=0 && i<=mxe);
|
---|
[5148] | 941 | i++,s0=1-(s1=directionge[i]),l0=l1;
|
---|
[3913] | 942 | }
|
---|
| 943 | on=ge[i];
|
---|
| 944 | if (i==tge)
|
---|
| 945 | s1=vg1;
|
---|
| 946 |
|
---|
| 947 | s =(ls-l0)/(l1-l0);
|
---|
| 948 | sg =s0*(1.0-s)+s*s1;
|
---|
| 949 | }
|
---|
[6412] | 950 | _assert_(on);
|
---|
[3913] | 951 | V.r= on->F(sg);
|
---|
| 952 | GV=VertexOnGeom(V,*on,sg);
|
---|
| 953 | return on;
|
---|
| 954 | }
|
---|
[12365] | 955 | /*}}}*/
|
---|
| 956 | /*FUNCTION Geometry::R2ToI2{{{*/
|
---|
[5180] | 957 | I2 Geometry::R2ToI2(const R2 & P) const {
|
---|
[5177] | 958 | /*coefIcoor is the coefficient used for integer coordinates:
|
---|
| 959 | * (x-pmin.x)
|
---|
| 960 | * Icoor x = (2^30 -1) ------------
|
---|
| 961 | * D
|
---|
| 962 | * where D is the longest side of the domain (direction x or y)
|
---|
| 963 | * so that (x-pmin.x)/D is in ]0 1[
|
---|
| 964 | *
|
---|
| 965 | * coefIcoor = (2^30 -1)/D
|
---|
| 966 | */
|
---|
| 967 | return I2( (Icoor1) (coefIcoor*(P.x-pmin.x)) ,(Icoor1) (coefIcoor*(P.y-pmin.y)) );
|
---|
[5130] | 968 | }/*}}}*/
|
---|
[12365] | 969 | /*FUNCTION Geometry::UnMarkEdges{{{*/
|
---|
[5130] | 970 | void Geometry::UnMarkEdges() {
|
---|
| 971 | for (int i=0;i<nbe;i++) edges[i].SetUnMark();
|
---|
| 972 | }/*}}}*/
|
---|
[3913] | 973 | }
|
---|