 |
Ice Sheet System Model
4.18
Code documentation
|
#include <Mesh.h>
|
| Mesh (BamgGeom *bamggeom, BamgMesh *bamgmesh, BamgOpts *bamgopts) |
|
| Mesh (int *index, double *x, double *y, int nods, int nels, BamgOpts *bamgopts) |
|
| Mesh (double *x, double *y, int nods, BamgOpts *bamgopts) |
|
| Mesh (Mesh &, Geometry *pGh=0, Mesh *pBTh=0, long maxnbv_in=0) |
|
| Mesh (const Mesh &, const int *flag, const int *bb, BamgOpts *bamgopts) |
|
| Mesh (long maxnbv, Mesh &BT, BamgOpts *bamgopts, int keepBackVertices=1) |
|
| Mesh (long maxnbv, Geometry &G, BamgOpts *bamgopts) |
|
| ~Mesh () |
|
const BamgVertex & | operator[] (long i) const |
|
BamgVertex & | operator[] (long i) |
|
const Triangle & | operator() (long i) const |
|
Triangle & | operator() (long i) |
|
void | SetIntCoor (const char *from=0) |
|
double | MinimalHmin () |
|
double | MaximalHmax () |
|
I2 | R2ToI2 (const R2 &P) const |
|
R2 | I2ToR2 (const I2 &P) const |
|
void | AddVertex (BamgVertex &s, Triangle *t, long long *=0) |
|
void | Insert (BamgOpts *bamgopts) |
|
void | Echo (void) |
|
void | ForceBoundary (BamgOpts *bamgopts) |
|
void | FindSubDomain (BamgOpts *bamgopts, int OutSide=0) |
|
long | TriangleReferenceList (long *) const |
|
void | TriangleIntNumbering (long *renumbering) |
|
void | CrackMesh (BamgOpts *bamgopts) |
|
void | SmoothMetric (BamgOpts *bamgopts, double raisonmax) |
|
void | BoundAnisotropy (BamgOpts *bamgopts, double anisomax, double hminaniso=1e-100) |
|
Edge ** | MakeGeomEdgeToEdge () |
|
long | SplitInternalEdgeWithBorderVertices () |
|
void | MakeBamgQuadtree () |
|
void | MaxSubDivision (BamgOpts *bamgopts, double maxsubdiv) |
|
void | NewPoints (Mesh &, BamgOpts *bamgopts, int KeepVertices=1) |
|
long | InsertNewPoints (long nbvold, long &NbTSwap, BamgOpts *bamgopts) |
|
void | TrianglesRenumberBySubDomain (bool justcompress=false) |
|
void | SmoothingVertex (BamgOpts *bamgopts, int=3, double=0.3) |
|
Metric | MetricAt (const R2 &) |
|
GeomEdge * | ProjectOnCurve (Edge &AB, BamgVertex &A, BamgVertex &B, double theta, BamgVertex &R, VertexOnEdge &BR, VertexOnGeom &GR) |
|
long | GetId (const Triangle &t) const |
|
long | GetId (const Triangle *t) const |
|
long | GetId (const BamgVertex &t) const |
|
long | GetId (const BamgVertex *t) const |
|
long | GetId (const Edge &t) const |
|
long | GetId (const Edge *t) const |
|
BamgVertex * | NearestVertex (int i, int j) |
|
Triangle * | TriangleFindFromCoord (const I2 &, long long[3], Triangle *tstart=0) |
|
void | ReadMesh (int *index, double *x, double *y, int nods, int nels, BamgOpts *bamgopts) |
|
void | ReadMesh (BamgMesh *bamgmesh, BamgOpts *bamgopts) |
|
void | WriteMesh (BamgMesh *bamgmesh, BamgOpts *bamgopts) |
|
void | ReadMetric (const BamgOpts *bamgopts) |
|
void | WriteMetric (BamgOpts *bamgopts) |
|
void | WriteIndex (int **pindex, int *pnels) |
|
void | AddMetric (BamgOpts *bamgopts) |
|
void | BuildMetric0 (BamgOpts *bamgopts) |
|
void | BuildMetric1 (BamgOpts *bamgopts) |
|
void | BuildGeometryFromMesh (BamgOpts *bamgopts=NULL) |
|
long | RandomNumber (long max) |
|
void | ReconstructExistingMesh (BamgOpts *bamgopts) |
|
void | CreateSingleVertexToTriangleConnectivity () |
|
void | UnMarkUnSwapTriangle () |
|
void | SetVertexFieldOn () |
|
void | SetVertexFieldOnBTh () |
|
|
void | TriangulateFromGeom1 (BamgOpts *bamgopts, int KeepVertices=1) |
|
void | TriangulateFromGeom0 (BamgOpts *bamgopts) |
|
void | Triangulate (double *x, double *y, int nods, BamgOpts *bamgopts) |
|
void | Init (long) |
|
int | ForceEdge (BamgVertex &a, BamgVertex &b, AdjacentTriangle &taret) |
|
int | SwapForForcingEdge (BamgVertex *&pva, BamgVertex *&pvb, AdjacentTriangle &tt1, long long &dets1, long long &detsa, long long &detsb, int &nbswap) |
|
Definition at line 21 of file Mesh.h.
◆ Mesh() [1/7]
Definition at line 13 of file Mesh.cpp.
13 :
Gh(*(
new Geometry())),
BTh(*
this){
28 if(bamggeom->
Edges==NULL) {
30 _printf_(
"WARNING: mesh present but no geometry found. Reconstructing...\n");
◆ Mesh() [2/7]
bamg::Mesh::Mesh |
( |
int * |
index, |
|
|
double * |
x, |
|
|
double * |
y, |
|
|
int |
nods, |
|
|
int |
nels, |
|
|
BamgOpts * |
bamgopts |
|
) |
| |
Definition at line 42 of file Mesh.cpp.
42 :
Gh(*(
new Geometry())),
BTh(*
this){
45 ReadMesh(index,x,y,nods,nels,bamgopts);
◆ Mesh() [3/7]
bamg::Mesh::Mesh |
( |
double * |
x, |
|
|
double * |
y, |
|
|
int |
nods, |
|
|
BamgOpts * |
bamgopts |
|
) |
| |
Definition at line 50 of file Mesh.cpp.
50 :
Gh(*(
new Geometry())),
BTh(*
this){
◆ Mesh() [4/7]
bamg::Mesh::Mesh |
( |
Mesh & |
Th, |
|
|
Geometry * |
pGh = 0 , |
|
|
Mesh * |
pBTh = 0 , |
|
|
long |
maxnbv_in = 0 |
|
) |
| |
Definition at line 153 of file Mesh.cpp.
154 :
Gh(*(pGh?pGh:&Th.Gh)),
BTh(*(pBth?pBth:
this)) {
157 maxnbv_in =
Max(maxnbv_in,Th.nbv);
199 triangles[i].Set(Th.triangles[i],Th,*
this);
201 edges[i].Set(Th,i,*
this);
203 vertices[i].Set(Th.vertices[i],Th,*
this);
◆ Mesh() [5/7]
bamg::Mesh::Mesh |
( |
const Mesh & |
Tho, |
|
|
const int * |
flag, |
|
|
const int * |
bb, |
|
|
BamgOpts * |
bamgopts |
|
) |
| |
Definition at line 54 of file Mesh.cpp.
54 :
Gh(*(
new Geometry())),
BTh(*
this) {
59 int * kk =
new int [Tho.nbv];
60 long * reft =
new long[Tho.nbt];
62 long * refv =
new long[Tho.nbv];
64 for (i=0;i<Tho.nbv;i++)
66 for (i=0;i<Tho.nbv;i++)
70 for (i=0;i<Tho.nbt;i++)
71 if( reft[i] >=0 && flag[i])
73 const Triangle & t = Tho.triangles[i];
75 kk[Tho.GetId(t[0])]=1;
76 kk[Tho.GetId(t[1])]=1;
77 kk[Tho.GetId(t[2])]=1;
78 itadj=Tho.GetId(t.TriangleAdj(0));
79 if ( reft[itadj] >=0 && !flag[itadj])
84 itadj=Tho.GetId(t.TriangleAdj(1));
85 if ( reft[itadj] >=0 && !flag[itadj])
89 itadj=Tho.GetId(t.TriangleAdj(2));
90 if ( reft[itadj] >=0 && !flag[itadj])
96 for (i=0;i<Tho.nbv;i++){
97 if (kk[i]>=0) kk[i]=k++;
99 _printf_(
" number of vertices " << k <<
", remove = " << Tho.nbv - k <<
"\n");
100 _printf_(
" number of triangles " << kt <<
", remove = " << nbInT-kt <<
"\n");
101 _printf_(
" number of New boundary edge " << nbNewBedge <<
"\n");
104 for (i=0;i<Tho.nbv;i++)
117 for (i=0;i<Tho.nbt;i++)
118 if(reft[i] >=0 && flag[i]){
119 const Triangle & t = Tho.triangles[i];
120 int i0 = Tho.GetId(t[0]);
121 int i1 = Tho.GetId(t[1]);
122 int i2 = Tho.GetId(t[2]);
123 if(i0<0 || i1<0 || i2<0){
125 _error_(
"i0<0 || i1<0 || i2< 0");
127 if(i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv){
129 _error_(
"i0>=Tho.nbv || i1>=Tho.nbv || i2>=Tho.nbv");
137 _error_(
"All triangles have been removed");
◆ Mesh() [6/7]
bamg::Mesh::Mesh |
( |
long |
maxnbv, |
|
|
Mesh & |
BT, |
|
|
BamgOpts * |
bamgopts, |
|
|
int |
keepBackVertices = 1 |
|
) |
| |
◆ Mesh() [7/7]
◆ ~Mesh()
◆ operator[]() [1/2]
const BamgVertex& bamg::Mesh::operator[] |
( |
long |
i | ) |
const |
|
inline |
◆ operator[]() [2/2]
◆ operator()() [1/2]
const Triangle& bamg::Mesh::operator() |
( |
long |
i | ) |
const |
|
inline |
◆ operator()() [2/2]
Triangle& bamg::Mesh::operator() |
( |
long |
i | ) |
|
|
inline |
◆ SetIntCoor()
void bamg::Mesh::SetIntCoor |
( |
const char * |
from = 0 | ) |
|
Definition at line 3678 of file Mesh.cpp.
3687 for (i=0;i<
nbv;i++) {
3698 long MaxICoord = 1073741823;
3701 _error_(
"coefIcoor should be positive, a problem in the geometry is likely");
3705 for (i=0;i<
nbv;i++) {
3710 int number_of_errors=0;
3711 for (i=0;i<
nbt;i++) {
3717 if (v0 && v1 && v2){
3720 triangles[i].
det=
det(v0->GetIntegerCoordinates(),v1->GetIntegerCoordinates(),v2->GetIntegerCoordinates());
3727 if (number_of_errors<20){
3737 if (number_of_errors)
_error_(
"Fatal error: some triangles have negative areas, see above");
◆ MinimalHmin()
double bamg::Mesh::MinimalHmin |
( |
| ) |
|
◆ MaximalHmax()
double bamg::Mesh::MaximalHmax |
( |
| ) |
|
◆ R2ToI2()
I2 bamg::Mesh::R2ToI2 |
( |
const R2 & |
P | ) |
const |
◆ I2ToR2()
R2 bamg::Mesh::I2ToR2 |
( |
const I2 & |
P | ) |
const |
◆ AddVertex()
Definition at line 1020 of file Mesh.cpp.
1039 long long det3local[3];
1045 BamgVertex* s0=t->GetVertex(0);
1046 BamgVertex* s1=t->GetVertex(1);
1047 BamgVertex* s2=t->GetVertex(2);
1050 long long detOld=t->det;
1055 int infvertexindex = s0 ? ((s1? (s2?-1:2):1)):0;
1058 if(((infvertexindex <0 ) && (detOld <0)) || ((infvertexindex >=0) && (detOld >0)) ){
1059 _error_(
"inconsistent configuration (Contact ISSM developers)");
1067 if (infvertexindex<0 ) {
1068 det3[0]=
bamg::det(s .GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s2->GetIntegerCoordinates());
1069 det3[1]=
bamg::det(s0->GetIntegerCoordinates(),s .GetIntegerCoordinates(),s2->GetIntegerCoordinates());
1070 det3[2]=
bamg::det(s0->GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s.GetIntegerCoordinates());}
1073 det3[0]= s0 ? -1 :
bamg::det(s.GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s2->GetIntegerCoordinates()) ;
1074 det3[1]= s1 ? -1 :
bamg::det(s0->GetIntegerCoordinates(),s.GetIntegerCoordinates(),s2->GetIntegerCoordinates()) ;
1075 det3[2]= s2 ? -1 :
bamg::det(s0->GetIntegerCoordinates(),s1->GetIntegerCoordinates(),s.GetIntegerCoordinates()) ;
1079 if (!det3[0]) izerodet=0,nbzerodet++;
1080 if (!det3[1]) izerodet=1,nbzerodet++;
1081 if (!det3[2]) izerodet=2,nbzerodet++;
1088 AdjacentTriangle ta = t->Adj(iedge);
1093 if (((Triangle*)ta)->
det<0 ) {
1101 _error_(
"Cannot add a vertex more than once. Check duplicates");
1106 t->SetUnMarkUnSwap(0);
1107 t->SetUnMarkUnSwap(1);
1108 t->SetUnMarkUnSwap(2);
1129 tt[0]->SetAdjAdj(0);
1130 tt[1]->SetAdjAdj(1);
1131 tt[2]->SetAdjAdj(2);
1137 tt[i0]->SetAdj2(i2,tt[i2],i0);
1138 tt[i1]->SetAdj2(i0,tt[i0],i1);
1139 tt[i2]->SetAdj2(i1,tt[i1],i2);
1141 tt[0]->SetSingleVertexToTriangleConnectivity();
1142 tt[1]->SetSingleVertexToTriangleConnectivity();
1143 tt[2]->SetSingleVertexToTriangleConnectivity();
1147 int rswap=tt[izerodet]->swap(iedge);
1150 _error_(
"swap the point s is on a edge");
◆ Insert()
void bamg::Mesh::Insert |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 2679 of file Mesh.cpp.
2689 if(bamgopts) verbose=bamgopts->
verbose;
2692 if (verbose>2)
_printf_(
" Insert initial " <<
nbv <<
" vertices\n");
2728 if (verbose>4)
_printf_(
" Prime Number = "<<PrimeNumber<<
"\n");
2729 if (verbose>4)
_printf_(
" k0 = "<<k0<<
"\n");
2732 for (i=0; i<
nbv; i++){
2742 if (++i>=
nbv)
_error_(
"all the vertices are aligned");
2744 if (verbose>4)
_printf_(
" i = "<<i<<
"\n");
2786 if (verbose>3)
_printf_(
" Begining of insertion process...\n");
2789 for (
long icount=2; icount<
nbv; icount++) {
2805 NbSwap += newvertex->Optim(1,0);
2810 _printf_(
" NbSwap of insertion: " << NbSwap <<
"\n");
◆ Echo()
void bamg::Mesh::Echo |
( |
void |
| ) |
|
Definition at line 2296 of file Mesh.cpp.
2305 for (i=0;i<
nbt;i++){
2306 _printf_(
" " << setw(4) << i+1 <<
": ["
2312 for (i=0;i<
nbv;i++){
◆ ForceBoundary()
void bamg::Mesh::ForceBoundary |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 2318 of file Mesh.cpp.
2323 int nbfe=0,nbswp=0,Nbswap=0;
2326 if(bamgopts) verbose=bamgopts->
verbose;
2329 if (verbose > 2)
_printf_(
" ForceBoundary nb of edge: " <<
nbe <<
"\n");
2332 for (
int t = 0; t <
nbt; t++){
2336 _error_(
"there is " << k <<
" triangles of mes = 0");
2340 AdjacentTriangle ta(0,0);
2341 for (
int i = 0; i <
nbe; i++){
2346 else Nbswap += nbswp;
2349 if ( nbswp < 0 && k < 5){
2355 _error_(
"There are " << k <<
" lost edges, the boundary might be crossing");
2357 for (
int j=0;j<
nbv;j++){
2360 if (verbose > 3)
_printf_(
" number of inforced edge = " << nbfe <<
", number of swap= " << Nbswap <<
"\n");
◆ FindSubDomain()
void bamg::Mesh::FindSubDomain |
( |
BamgOpts * |
bamgopts, |
|
|
int |
OutSide = 0 |
|
) |
| |
Definition at line 2363 of file Mesh.cpp.
2369 if(bamgopts) verbose=bamgopts->
verbose;
2372 if (OutSide)
_printf_(
" Find all external sub-domain\n");
2373 else _printf_(
" Find all internal sub-domain\n");
2375 short * HeapArete =
new short[
nbt];
2376 Triangle ** HeapTriangle =
new Triangle* [
nbt];
2383 long NbSubDomTot =0;
2384 for(it=0;it<
nbt;it++) {
2391 HeapTriangle[i] =t ;
2395 {
while (HeapArete[i]--)
2397 int na = HeapArete[i];
2398 Triangle * tc = HeapTriangle[i];
2399 if( ! tc->Locked(na))
2408 HeapTriangle[i] = ta;
2438 delete [] HeapArete;
2439 delete [] HeapTriangle;
2440 _error_(
"The boundary is not close: all triangles are outside");
2443 delete [] HeapArete;
2444 delete [] HeapTriangle;
2456 long * mark =
new long[
nbt];
2457 for (it=0;it<
nbt;it++)
2463 if (mark[it] == -1) {
2485 for (it=0;it<
nbt && nbk ;it++)
2486 for (
int na=0;na<3 && nbk ;na++)
2489 long kl = ta ? mark[
GetId(ta)] : -2;
2492 if (kl >=0 &&
subdomains[kl].ReferenceNumber <0 && kr >=0 &&
subdomains[kr].ReferenceNumber>=0)
2494 if (kr >=0 &&
subdomains[kr].ReferenceNumber <0 && kl >=0 &&
subdomains[kl].ReferenceNumber>=0)
2496 if(kr<0 && kl >=0 &&
subdomains[kl].ReferenceNumber>=0)
2498 if(kl<0 && kr >=0 &&
subdomains[kr].ReferenceNumber>=0)
2517 if(verbose>4)
_printf_(
" Number of removes subdomains (OutSideMesh) = " <<
nbsubdomains-j <<
"\n");
2531 long * mark =
new long[
nbt];
2534 for (it=0;it<
nbt;it++)
2542 Edge &e = *GeomEdgetoEdge[
Gh.
GetId(eg)];
2544 BamgVertex * v0 = e(0),*v1 = e(1);
2545 Triangle *t = v0->t;
2548 if (((eg[0].r-eg[1].r),(e[0].r-e[1].r))<0) direction = -direction ;
2557 if (ta.EdgeVertex(0) == v1) {
2562 if(t<triangles || t >=
triangles+
nbt || t->det < 0 || t->link == 0) {
2563 _error_(
"bad definition of SubSomain " << i);
2577 if (mark[
GetId(tt)]>=0){
2578 _error_(
"mark[GetId(tt)]>=0");
2586 if(t == (Triangle *) ta) {
2587 _error_(
"bad definition of SubSomain " << i);
2593 if (verbose>5)
_printf_(
"WARNING: " <<
nbsubdomains-inew <<
" SubDomains are being removed\n");
2596 for (it=0;it<
nbt;it++)
2597 if ( mark[it] ==-1 )
2599 delete [] GeomEdgetoEdge;
2604 for (it=0;it<
nbt;it++)
◆ TriangleReferenceList()
long bamg::Mesh::TriangleReferenceList |
( |
long * |
reft | ) |
const |
Definition at line 4046 of file Mesh.cpp.
4053 for (
int it=0;it<
nbt;it++) reft[it]=-1;
4062 if (!t0){
_error_(
"At least one subdomain is empty");}
4077 }
while (t0 != (t=t->link));
◆ TriangleIntNumbering()
void bamg::Mesh::TriangleIntNumbering |
( |
long * |
renumbering | ) |
|
Definition at line 4036 of file Mesh.cpp.
4039 for (
int i=0;i<
nbt;i++){
4041 else renumbering[i]=-1;
◆ CrackMesh()
void bamg::Mesh::CrackMesh |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 2148 of file Mesh.cpp.
2152 int i,j,k,num,count;
2158 if(bamgopts) verbose=bamgopts->
verbose;
2161 for (k=i=0;i<
nbe;i++){
2162 if(
edges[i].GeomEdgeHook->Cracked()) k++;
2167 if (verbose>4)
_printf_(
" number of Cracked Edges = " << k <<
"\n");
2177 int* splitvertex=
new int[
nbv];
2178 for (i=0;i<
nbv;i++) splitvertex[i]=0;
2180 for (i=0;i<
nbe;i++){
2181 if(
edges[i].GeomEdgeHook->Cracked()){
2199 if (splitvertex[i1]==3 || splitvertex[i2]==3){
2200 delete [] splitvertex;
2201 _error_(
"Crossing rifts not supported yet");
2212 for (i=0;i<
nbv;i++){
2213 if (splitvertex[i]==2){
2221 delete [] splitvertex;
2246 if (
GetId((*ta.t)[j])==i1){
2253 if (
GetId((*ta.t)[j])==i2){
2261 if (Edgeflags[i]==0){
2279 if (count++>50)
_error_(
"Maximum number of iteration exceeded");
2280 }
while ((tbegin != ta));
2285 if (Edgeflags[i]!=2){
2286 _error_(
"A problem occured: at least one crack edge (number " << i+1 <<
") does not belong to 2 elements");
2289 delete [] Edgeflags;
◆ SmoothMetric()
void bamg::Mesh::SmoothMetric |
( |
BamgOpts * |
bamgopts, |
|
|
double |
raisonmax |
|
) |
| |
Definition at line 3785 of file Mesh.cpp.
3792 if(bamgopts) verbose=bamgopts->
verbose;
3794 if(raisonmax<1.1)
return;
3795 if(verbose > 1)
_printf_(
" Mesh::SmoothMetric raisonmax = " << raisonmax <<
"\n");
3798 long *first_np_or_next_t0 =
new long[
nbv];
3799 long *first_np_or_next_t1 =
new long[
nbv];
3800 long Head0 =0,Head1=-1;
3801 double logseuil= log(raisonmax);
3803 for(i=0;i<
nbv-1;i++)
3804 first_np_or_next_t0[i]=i+1;
3805 first_np_or_next_t0[
nbv-1]=-1;
3807 first_np_or_next_t1[i]=-1;
3809 while(Head0>=0&& kk++<100){
3811 for(i=Head0;i>=0;i=first_np_or_next_t0[ip=i],first_np_or_next_t0[ip]=-1) {
3819 BamgVertex *pvj0 = ta.EdgeVertex(0);
3822 if (
vertices+i != ta.EdgeVertex(1)){
3823 _error_(
"vertices+i != ta.EdgeVertex(1)");
3825 BamgVertex *pvj = (ta.EdgeVertex(0));
3826 BamgVertex & vj = *pvj;
3829 if (j<0 || j >=
nbv){
3832 R2 Aij = (
R2) vj - (
R2) vi;
3835 double hi = ll/vi.m.Length(Aij.x,Aij.y);
3836 double hj = ll/vj.m.Length(Aij.x,Aij.y);
3839 double dh=(hj-hi)/ll;
3841 vj.m.IntersectWith(vi.m/(1 +logseuil*ll/hi));
3842 if(first_np_or_next_t1[j]<0)
3843 kch++,first_np_or_next_t1[j]=Head1,Head1=j;
3849 double li = vi.m.Length(Aij.x,Aij.y);
3850 if( vj.m.IntersectWith(vi.m/(1 +logseuil*li)) )
3851 if(first_np_or_next_t1[j]<0)
3852 kch++,first_np_or_next_t1[j]=Head1,Head1=j;
3855 if ( &vj == pvj0 )
break;
3860 Exchange(first_np_or_next_t0,first_np_or_next_t1);
3862 if(verbose>2)
_printf_(
" number of iterations = " << kch <<
"\n");
3863 delete [] first_np_or_next_t0;
3864 delete [] first_np_or_next_t1;
◆ BoundAnisotropy()
void bamg::Mesh::BoundAnisotropy |
( |
BamgOpts * |
bamgopts, |
|
|
double |
anisomax, |
|
|
double |
hminaniso = 1e-100 |
|
) |
| |
Definition at line 1155 of file Mesh.cpp.
1159 double lminaniso = 1./ (
Max(hminaniso*hminaniso,1e-100));
1162 if(bamgopts) verbose=bamgopts->
verbose;
1165 if (verbose > 1)
_printf_(
" BoundAnisotropy by " << anisomax <<
"\n");
1167 double h1=1.e30,h2=1e-30;
1168 double coef = 1./(anisomax*anisomax);
1169 double hn1=1.e30,hn2=1e-30,rnx =1.e-30,rx=0;
1172 for (
int i=0;i<
nbv;i++){
1174 double lmax=Vp.lmax();
1175 Vp*=
Min(lminaniso,lmax)/lmax;
1176 Vp.BoundAniso2(coef);
1181 h1 =
Min(h1,Vp.lmin());
1182 h2 =
Max(h2,Vp.lmax());
1183 hn1=
Min(hn1,Vp.lmin());
1184 hn2=
Max(hn2,Vp.lmax());
1185 rx =
Max(rx,Vp.Aniso2());
1186 rnx=
Max(rnx,Vp.Aniso2());
1192 _printf_(
" input: Hmin = " << pow(h2,-0.5) <<
", Hmax = " << pow(h1,-0.5) <<
", factor of anisotropy max = " << pow(rx,0.5) <<
"\n");
1193 _printf_(
" output: Hmin = " << pow(hn2,-0.5) <<
", Hmax = " << pow(hn1,-0.5)<<
", factor of anisotropy max = " <<pow(rnx,0.5) <<
"\n");
◆ MakeGeomEdgeToEdge()
Edge ** bamg::Mesh::MakeGeomEdgeToEdge |
( |
| ) |
|
Definition at line 2899 of file Mesh.cpp.
2905 Edge **e=
new Edge* [
Gh.
nbe];
2908 for ( i=0;i<
Gh.
nbe ; i++)
2910 for ( i=0;i<
nbe ; i++)
2912 Edge * ei =
edges+i;
2914 e[
Gh.
GetId(GeomEdgeHook)] = ei;
2916 for ( i=0;i<
nbe ; i++)
2917 for (
int ii=0;ii<2;ii++) {
2918 Edge * ei =
edges+i;
2921 while (!(*GeomEdgeHook)[j].Required()) {
2922 Adj(GeomEdgeHook,j);
2924 if (e[
Gh.
GetId(GeomEdgeHook)])
break;
2925 e[
Gh.
GetId(GeomEdgeHook)] = ei;
2930 for ( i=0;i<
Gh.
nbe ; i++){
2933 if(kk<10)
_printf_(
"BUG: the geometrical edge " << i <<
" is on no edge curve\n");
◆ SplitInternalEdgeWithBorderVertices()
long bamg::Mesh::SplitInternalEdgeWithBorderVertices |
( |
| ) |
|
Definition at line 3867 of file Mesh.cpp.
3875 for (it=0;it<
nbt;it++){
3878 for (
int j=0;j<3;j++)
3879 if(!t.Locked(j) && !t.Hidden(j)){
3882 Triangle &tt = *ptt;
3884 if( ptt && tt.link && it <
GetId(tt))
3888 if (v0.GeomEdgeHook && v1.GeomEdgeHook){
3889 R2 P= ((
R2) v0 + (
R2) v1)*0.5;
3905 for (
int i=nbvold;i<
nbv;i++) {
3911 vi.ReferenceNumber=0;
3913 if (tcvi && !tcvi->link) {
3914 _printf_(
"problem inserting point in SplitInternalEdgeWithBorderVertices (tcvj && !tcvj->link)\n");
3918 if (!tcvi || tcvi->det<0){
3919 _error_(
"!tcvi || tcvi->det < 0");
3922 NbSwap += vi.Optim(1);
3926 _printf_(
" number of points: " << iv <<
"\n");
3927 _printf_(
" number of swap to split internal edges with border vertices: " << NbSwap <<
"\n");
3931 if (NbSplitEdge>
nbv-nbvold)
_printf_(
"WARNING: not enough vertices to split all internal edges, we lost " << NbSplitEdge - (
nbv-nbvold) <<
" edges...\n");
3932 if (verbose>2)
_printf_(
"SplitInternalEdgeWithBorderVertices: Number of splited edge " << NbSplitEdge <<
"\n");
◆ MakeBamgQuadtree()
void bamg::Mesh::MakeBamgQuadtree |
( |
| ) |
|
◆ MaxSubDivision()
void bamg::Mesh::MaxSubDivision |
( |
BamgOpts * |
bamgopts, |
|
|
double |
maxsubdiv |
|
) |
| |
Definition at line 2953 of file Mesh.cpp.
2960 if(bamgopts) verbose=bamgopts->
verbose;
2962 const double maxsubdiv2 = maxsubdiv*maxsubdiv;
2963 if(verbose>1)
_printf_(
" Limit the subdivision of a edges in the new mesh by " << maxsubdiv <<
"\n");
2968 for (it=0;it<
nbt;it++){
2970 for (
int j=0;j<3;j++){
2972 Triangle &tt = *ptt;
2973 if ( (!ptt || it <
GetId(tt)) && ( tt.link || t.link)){
2978 double l = M(AB,AB);
2981 R2 AC = M.Orthogonal(AB);
2982 double lc = M(AC,AC);
2984 D2xD2 Rt1(Rt.inv());
2985 D2xD2 D(maxsubdiv2,0,0,lc);
2986 D2xD2 MM = Rt1*D*Rt1.t();
2987 v0.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);
2994 R2 AC = M.Orthogonal(AB);
2995 double lc = M(AC,AC);
2997 D2xD2 Rt1(Rt.inv());
2998 D2xD2 D(maxsubdiv2,0,0,lc);
2999 D2xD2 MM = Rt1*D*Rt1.t();
3000 v1.m = M = Metric(MM.x.x,MM.y.x,MM.y.y);
3007 _printf_(
" number of metric changes = " << nbchange <<
", maximum number of subdivision of a edges before change = " << pow(lmax,0.5) <<
"\n");
◆ NewPoints()
void bamg::Mesh::NewPoints |
( |
Mesh & |
Bh, |
|
|
BamgOpts * |
bamgopts, |
|
|
int |
KeepVertices = 1 |
|
) |
| |
Definition at line 3040 of file Mesh.cpp.
3050 long* first_np_or_next_t=
new long[
maxnbt];
3054 if(bamgopts) verbose=bamgopts->
verbose;
3057 if(KeepVertices && (&Bh !=
this) && (
nbv+Bh.nbv<
maxnbv)){
3058 if (verbose>5)
_printf_(
" Inserting initial mesh points\n");
3059 bool pointsoutside =
false;
3060 for(i=0;i<Bh.nbv;i++){
3061 BamgVertex &bv=Bh[i];
3065 if(tcvj->det<0 || !tcvj->link){
3066 pointsoutside =
true;
3069 IssmDouble area_1=((bv.r.x -(*tcvj)(2)->r.x)*((*tcvj)(1)->r.y-(*tcvj)(2)->r.y)
3070 - (bv.r.y -(*tcvj)(2)->r.y)*((*tcvj)(1)->r.x-(*tcvj)(2)->r.x));
3071 IssmDouble area_2=(((*tcvj)(0)->r.x -(*tcvj)(2)->r.x)*(bv.r.y -(*tcvj)(2)->r.y)
3072 - ((*tcvj)(0)->r.y -(*tcvj)(2)->r.y)*(bv.r.x -(*tcvj)(2)->r.x));
3073 IssmDouble area_3 =((bv.r.x -(*tcvj)(1)->r.x)*((*tcvj)(0)->r.y-(*tcvj)(1)->r.y)
3074 - (bv.r.y -(*tcvj)(1)->r.y)*((*tcvj)(0)->r.x-(*tcvj)(1)->r.x));
3075 if(area_1<0 || area_2<0 || area_3<0){
3076 pointsoutside =
true;
3079 if(!bv.GeomEdgeHook){
3086 Bh.CreateSingleVertexToTriangleConnectivity();
3089 else Bh.CreateSingleVertexToTriangleConnectivity();
3092 for(i=0;i<
nbt;i++) first_np_or_next_t[i]=-(i+1);
3097 if (verbose>5)
_printf_(
" Big loop\n");
3106 next_t=-first_np_or_next_t[i];
3110 if (i<0 || i>=
nbt ){
3111 _error_(
"Index problem in NewPoints (i=" << i <<
" not in [0 " <<
nbt-1 <<
"])");
3114 first_np_or_next_t[i] = iter;
3118 AdjacentTriangle tj(t,j);
3119 BamgVertex &vA = *tj.EdgeVertex(0);
3120 BamgVertex &vB = *tj.EdgeVertex(1);
3123 if (!t->link)
continue;
3124 if (t->det <0)
continue;
3125 if (t->Locked(j))
continue;
3127 AdjacentTriangle tadjj = t->Adj(j);
3131 if (ta->det<0)
continue;
3138 if(first_np_or_next_t[k]==iter)
continue;
3146 for (i=nbtold;i<
nbt;i++) first_np_or_next_t[i]=iter;
3150 for (i=nbvold;i<
nbv;i++){
3153 Triangle* tbegin= (Triangle*) ta;
3156 kt =
GetId((Triangle*) ta);
3157 if (first_np_or_next_t[kt]>0){
3158 first_np_or_next_t[kt]=-Headt;
3161 if (ta.EdgeVertex(0)!=s){
3162 _error_(
"ta.EdgeVertex(0)!=s");
3165 }
while ( (tbegin != (Triangle*) ta));
3168 }
while(
nbv!=nbvold);
3169 delete [] first_np_or_next_t;
3172 for(i=0;i<
nbv;i++) NbSwapf +=
vertices[i].Optim(0);
◆ InsertNewPoints()
long bamg::Mesh::InsertNewPoints |
( |
long |
nbvold, |
|
|
long & |
NbTSwap, |
|
|
BamgOpts * |
bamgopts |
|
) |
| |
Definition at line 2815 of file Mesh.cpp.
2819 double seuil= 1.414/2.;
2825 if(bamgopts) verbose=bamgopts->
verbose;
2828 const long nbvnew=
nbv-nbvold;
2831 if (verbose>5)
_printf_(
" Try to Insert " << nbvnew <<
" new points\n");
2834 if (!nbvnew)
return 0;
2840 for (
int is3=0; is3<nbvnew; is3++){
2841 long j=nbvold +(k3 = (k3+PrimeNumber)%nbvnew);
2849 for(i=nbvold;i<
nbv;i++){
2861 _error_(
"&vj!= orderedvertices[j]");
2867 vj.ReferenceNumber=0;
2869 if (tcvj && !tcvj->link){
2870 _printf_(
"While trying to add the following point:\n");
2872 _printf_(
"BAMG determined that it was inside the following triangle, which probably lies outside of the geometric domain\n");
2874 _error_(
"problem inserting point in InsertNewPoints (tcvj=" << tcvj <<
" and tcvj->link=" << tcvj->link <<
")");
2878 NbSwap += vj.Optim(1);
2882 vi.PreviousNumber = 0;
2886 _printf_(
" number of new points: " << iv <<
"\n");
2887 _printf_(
" number of to close (?) points: " <<
nbv-iv <<
"\n");
2888 _printf_(
" number of swap after: " << NbSwap <<
"\n");
2892 for (i=nbvold;i<
nbv;i++) NbSwap +=
vertices[i].Optim(1);
2893 if (verbose>3)
_printf_(
" NbSwap=" << NbSwap <<
"\n");
◆ TrianglesRenumberBySubDomain()
void bamg::Mesh::TrianglesRenumberBySubDomain |
( |
bool |
justcompress = false | ) |
|
Definition at line 3611 of file Mesh.cpp.
3614 long *renu=
new long[
nbt];
3618 for ( it=0;it<
nbt;it++)
3628 if (kt<0 || kt >=
nbt ){
3636 while (t0 != (t=t->link));
3640 for ( k=0,it=0;it<
nbt;it++)
3645 for ( it=0;it<
nbt;it++){
3646 if (renu[it]==-1) renu[it]=k++;
3652 for ( it=0;it<
nbt;it++)
3660 for ( it=0;it<
nbt;it++)
3665 while ( (j=renu[i]) >= 0)
◆ SmoothingVertex()
void bamg::Mesh::SmoothingVertex |
( |
BamgOpts * |
bamgopts, |
|
|
int |
nbiter = 3 , |
|
|
double |
omega = 0.3 |
|
) |
| |
Definition at line 3740 of file Mesh.cpp.
3747 if(bamgopts) verbose=bamgopts->
verbose;
3754 Triangle ** tstart=
new Triangle* [
nbv];
3758 for ( i=0;i<
nbv;i++)
3761 for ( i=0;i<
nbv;i++)
3767 if(verbose>2)
_printf_(
" SmoothingVertex: nb Iteration = " << nbiter <<
", Omega=" << omega <<
"\n");
3768 for (k=0;k<nbiter;k++)
3772 for ( i=0;i<
nbv;i++)
3773 if (tstart[i] != &vide)
3774 delta=
Max(delta,
vertices[i].Smoothing(*
this,
BTh,tstart[i],omega));
3775 for ( i=0;i<
nbv;i++)
3776 if (tstart[i] != &vide)
3778 if (verbose>3)
_printf_(
" move max = " << pow(delta,0.5) <<
", iteration = " << k <<
", nb of swap = " << NbSwap <<
"\n");
◆ MetricAt()
Metric bamg::Mesh::MetricAt |
( |
const R2 & |
A | ) |
|
Definition at line 3011 of file Mesh.cpp.
3020 return Metric(ba,*edge.EdgeVertex(0),bb,*edge.EdgeVertex(1));}
3023 double s = deta[0]+deta[1]+deta[2];
3027 return Metric(aa,(*t)[0],(*t)[1],(*t)[2]);
◆ ProjectOnCurve()
Definition at line 3174 of file Mesh.cpp.
3181 BamgVertex * pvA=&vA, * pvB=&vB;
3183 pA=vA.BackgroundVertexHook;
3186 pA=vA.BackgroundEdgeHook->be;
3187 tA=vA.BackgroundEdgeHook->abcisse;
3190 _error_(
"ProjectOnCurve On BamgVertex " <<
BTh.
GetId(vA) <<
" forget call to SetVertexFieldOnBTh");
3194 pB=vB.BackgroundVertexHook;
3197 pB=vB.BackgroundEdgeHook->be;
3198 tB=vB.BackgroundEdgeHook->abcisse;
3201 _error_(
"ProjectOnCurve On BamgVertex " <<
BTh.
GetId(vB) <<
" forget call to SetVertexFieldOnBTh");
3204 if (!pA || !pB || !e){
3211 _error_(
"e<BTh.edges || e>=BTh.edges+BTh.nbe");
3220 e = vA.BackgroundEdgeHook->be;
3230 e = vB.BackgroundEdgeHook->be;
3234 _error_(
"case not supported yet");
3240 double cosE01AB = (( (
R2) (*e)[1] - (
R2) (*e)[0] ) , AB);
3242 int direction = (cosE01AB>0) ? 1 : 0;
3245 double abscisse = -1;
3247 for (
int step=0;step<2;step++){
3252 BamgVertex *v0=pvA,*v1;
3257 for ( eee=e,iii=direction,te0=tA;
3258 eee && (((
void*) eee) != pB) && ((
void*) (v1=&((*eee)[iii]))) != pB ;
3259 neee = eee->adj[iii],iii = 1-neee->Intersection(*eee),eee = neee,v0=v1,te0=1-iii ) {
3267 if (step && abscisse <= lg) {
3268 double sss = (abscisse-lg0)/dp;
3269 double thetab = te0*(1-sss)+ sss*iii;
3271 BR = VertexOnEdge(&
R,eee,thetab);
3277 if ((
void*) v1 == pB)
3285 abscisse = lg*theta;
3286 if (abscisse <= lg && abscisse >= lg0 )
3288 double sss = (abscisse-lg0)/dp;
3289 double thetab = te0*(1-sss)+ sss*tB;
3291 BR = VertexOnEdge(&
R,eee,thetab);
3295 abscisse = lg*theta;
◆ GetId() [1/6]
long bamg::Mesh::GetId |
( |
const Triangle & |
t | ) |
const |
◆ GetId() [2/6]
long bamg::Mesh::GetId |
( |
const Triangle * |
t | ) |
const |
◆ GetId() [3/6]
long bamg::Mesh::GetId |
( |
const BamgVertex & |
t | ) |
const |
◆ GetId() [4/6]
long bamg::Mesh::GetId |
( |
const BamgVertex * |
t | ) |
const |
◆ GetId() [5/6]
long bamg::Mesh::GetId |
( |
const Edge & |
t | ) |
const |
◆ GetId() [6/6]
long bamg::Mesh::GetId |
( |
const Edge * |
t | ) |
const |
◆ NearestVertex()
BamgVertex * bamg::Mesh::NearestVertex |
( |
int |
i, |
|
|
int |
j |
|
) |
| |
◆ TriangleFindFromCoord()
Triangle * bamg::Mesh::TriangleFindFromCoord |
( |
const I2 & |
B, |
|
|
long long |
det3[3], |
|
|
Triangle * |
tstart = 0 |
|
) |
| |
Definition at line 3945 of file Mesh.cpp.
3953 if (tstart) t=tstart;
3959 if (!
quadtree)
_error_(
"no starting triangle provided and no quadtree available");
3965 if (!a)
_error_(
"problem while trying to find nearest vertex from a given point. No output found");
3966 if (!a->t)
_error_(
"no triangle is associated to vertex number " <<
GetId(a)+1 <<
" (orphan?)");
3983 int k0=(*t)(0)?(((*t)(1)?((*t)(2)?-1:2):1)):0;
3986 det3[k0]=
det(B,(*t)[k1],(*t)[k2]);
3987 det3[k1]=det3[k2]=-1;
4001 if (++counter>=10000)
_error_(
"Maximum number of iteration reached (threshold = " << counter <<
").");
4007 det3[jp]=
det(*(*t)(j),*(*t)(jn),B);
4008 det3[jn] = t->det-det3[j] -det3[jp];
4012 if (det3[0] < 0 ) ii[k++]=0;
4013 if (det3[1] < 0 ) ii[k++]=1;
4014 if (det3[2] < 0 ) ii[k++]=2;
4022 AdjacentTriangle t1 = t->Adj(jj=ii[0]);
4023 if ((t1.det() < 0 ) && (k == 2))
4024 t1 = t->Adj(jj=ii[1]);
◆ ReadMesh() [1/2]
void bamg::Mesh::ReadMesh |
( |
int * |
index, |
|
|
double * |
x, |
|
|
double * |
y, |
|
|
int |
nods, |
|
|
int |
nels, |
|
|
BamgOpts * |
bamgopts |
|
) |
| |
Definition at line 249 of file Mesh.cpp.
255 bool* nodeflags=NULL;
260 if(bamgopts) verbose=bamgopts->
verbose;
263 if (verbose)
_printf_(
"Reading vertices (" <<
nbv <<
")\n");
276 if (verbose)
_printf_(
"Reading triangles (" <<
nbt <<
")\n");
278 nodeflags=xNew<bool>(
nbv);
279 for(i=0;i<
nbv;i++) nodeflags[i]=
false;
283 i1=(long)index[i*3+0]-1;
284 i2=(long)index[i*3+1]-1;
285 i3=(long)index[i*3+2]-1;
286 t=Triangle(
this,i1,i2,i3);
288 nodeflags[i1]=nodeflags[i2]=nodeflags[i3]=
true;
292 if (verbose)
_printf_(
"Building Geometry\n");
294 if (verbose)
_printf_(
"Completing geometry\n");
301 _printf_(
"Vertex " << i+1 <<
" does not belong to any element\n");
305 if(isorphan)
_error_(
"Orphan found in mesh, see ids above");
308 xDelete<bool>(nodeflags);
◆ ReadMesh() [2/2]
Definition at line 311 of file Mesh.cpp.
314 double Hmin = HUGE_VAL;
322 if(bamgopts) verbose=bamgopts->
verbose;
330 if(verbose>5)
_printf_(
" processing Vertices\n");
345 if(verbose>5)
_error_(
"no Vertices found in the initial mesh");
350 if(verbose>5)
_printf_(
" processing Triangles\n");
358 t=Triangle(
this,i1,i2,i3);
359 t.color=(long)bamgmesh->
Triangles[i*4+3];
363 if(verbose>5)
_error_(
"no Triangles found in the initial mesh");
368 if(verbose>5)
_printf_(
" processing VerticesOnGeomEdge\n");
383 if(verbose>5)
_printf_(
" processing VerticesOnGeomVertex\n");
395 if (bamgmesh->
Edges){
399 if(verbose>5)
_printf_(
" processing Edges\n");
403 len=
new double[
nbv];
404 for(i=0;i<
nbv;i++) len[i]=0;
407 i1=(int)bamgmesh->
Edges[i*3+0]-1;
408 i2=(
int)bamgmesh->
Edges[i*3+1]-1;
415 double l12=sqrt((x12,x12));
422 Hmin =
Min(Hmin,l12);
440 BamgVertex *v=
edges[i].
v[j];
459 if(verbose>5)
_printf_(
" processing EdgesOnGeomEdge\n");
462 for (i1=0;i1<i2;i1++) {
466 if(!(i>=0 && j>=0 && i<
nbe && j<
Gh.
nbe)) {
467 _error_(
"ReadMesh error: EdgesOnGeomEdge edge provided (line " << i1+1 <<
": [" << i+1 <<
" " << j+1 <<
"]) is incorrect (must be positive, [0<i<nbe=" <<
nbe <<
" 0<j<Gh.nbe=" <<
Gh.
nbe <<
"]");
475 long i3,head,direction;
476 if(verbose>5)
_printf_(
" processing SubDomains\n");
483 if (i3!=3)
_error_(
"Bad Subdomain definition: first number should be 3");
484 if (head<0 || head>=
nbt)
_error_(
"Bad Subdomain definition: head should in [1 " <<
nbt <<
"] (triangle number)");
◆ WriteMesh()
Definition at line 493 of file Mesh.cpp.
500 int* connectivitysize_1=NULL;
501 int connectivitymax_1=0;
508 if(bamgopts) verbose=bamgopts->
verbose;
511 long* reft =
new long[
nbt];
512 long* numt =
new long[
nbt];
519 head_1=xNew<int>(
nbv);
520 next_1=xNew<int>(3*
nbt);
521 connectivitysize_1=xNew<int>(
nbv);
524 for (i=0;i<
nbv;i++) head_1[i]=-1;
525 for (i=0;i<
nbv;i++) connectivitysize_1[i]=0;
528 for (i=0;i<
nbt;i++) {
533 if (k>3*
nbt-1 || k<0)
_error_(
"k = " << k <<
", nbt = " <<
nbt);
535 if (v>
nbv-1 || v<0)
_error_(
"v = " << v <<
", nbv = " <<
nbv);
537 connectivitysize_1[v]+=1;
544 if (connectivitysize_1[i]>connectivitymax_1) connectivitymax_1=connectivitysize_1[i];
550 if(verbose>5)
_printf_(
" writing Vertices\n");
565 if(verbose>5)
_printf_(
" writing Edges\n");
568 int NumIssmSegments=0;
570 bamgmesh->
Edges=xNew<double>(3*
nbe);
575 if(
edges[i].GeomEdgeHook){
582 if(verbose>5)
_printf_(
" writing element edges\n");
583 SetOfEdges4* edge4=
new SetOfEdges4(
nbt*3,
nbv);
584 double* elemedge=NULL;
585 elemedge=xNew<double>(3*
nbt);
586 for (i=0;i<3*
nbt;i++) elemedge[i]=-2.;
594 n =edge4->SortAndFind(i1,i2);
597 n=edge4->SortAndAdd(i1,i2);
598 elemedge[n*2+0]=double(k);
602 elemedge[n*2+1]=double(k);
610 bamgmesh->
IssmEdges=xNew<double>(4*edge4->nb());
611 for (i=0;i<edge4->nb();i++){
618 bamgmesh->
IssmEdges[i*4+0]=edge4->i(i)+1;
619 bamgmesh->
IssmEdges[i*4+1]=edge4->j(i)+1;
622 bamgmesh->
IssmEdges[i*4+0]=edge4->j(i)+1;
623 bamgmesh->
IssmEdges[i*4+1]=edge4->i(i)+1;
630 bamgmesh->
IssmEdges[i*4+2]=elemedge[2*i+0]+1;
631 bamgmesh->
IssmEdges[i*4+3]=elemedge[2*i+1]+1;
635 xDelete<double>(elemedge);
638 if(verbose>5)
_printf_(
" writing IssmSegments\n");
644 if(
edges[i].GeomEdgeHook){
649 for(j=head_1[i1];j!=-1;j=next_1[j]){
675 _error_(
"Element holding segment [" << i1+1 <<
" " << i2+1 <<
"] not found...");
681 if(verbose>5)
_printf_(
" writing Triangles\n");
691 if (reft[i]>=0 && !( t.Hidden(0) || t.Hidden(1) || t.Hidden(2) )){
702 if(verbose>5)
_printf_(
" writing SubDomains\n");
716 if(verbose>5)
_printf_(
" writing SubDomainsFromGeom\n");
730 if(verbose>5)
_printf_(
" writing VerticesOnGeomVertex\n");
744 if(verbose>5)
_printf_(
" writing VerticesOnGeomEdge\n");
751 if (!v.OnGeomEdge()){
752 _error_(
"A vertices supposed to be OnGeomEdge is actually not");
761 if(verbose>5)
_printf_(
" writing EdgesOnGeomEdge\n");
764 if (
edges[i].GeomEdgeHook) k=k+1;
772 if (
edges[i].GeomEdgeHook){
781 if(verbose>5)
_printf_(
" writing Element connectivity\n");
801 if(verbose>5)
_printf_(
" writing Nodal element connectivity\n");
808 for(j=head_1[i];j!=-1;j=next_1[j]){
809 _assert_(connectivitymax_1*i+k < connectivitymax_1*
nbv);
816 if(verbose>5)
_printf_(
" writing Nodal connectivity\n");
820 int* connectivitysize_2=NULL;
821 int connectivitymax_2=0;
824 head_2=xNew<int>(
nbv);
825 next_2=xNew<int>(2*i1);
826 connectivitysize_2=xNew<int>(
nbv);
828 for (i=0;i<
nbv;i++) head_2[i]=-1;
829 for (i=0;i<
nbv;i++) connectivitysize_2[i]=0;
834 int v=(int)bamgmesh->
IssmEdges[i*i2+j]-1;
835 if (k>2*i1-1 || k<0)
_error_(
"Index exceed matrix dimensions (k=" << k <<
" not in [0 " << 2*i1-1 <<
"]");
837 if (v>
nbv-1 || v<0)
_error_(
"Index exceed matrix dimensions (v=" << v <<
" not in [0 " <<
nbv-1 <<
"])");
839 connectivitysize_2[v]+=1;
844 if (connectivitysize_2[i]>connectivitymax_2) connectivitymax_2=connectivitysize_2[i];
854 for(j=head_2[i];j!=-1;j=next_2[j]){
855 _assert_(connectivitymax_2*i+k < connectivitymax_2*
nbv);
856 num=(int)bamgmesh->
IssmEdges[
int(j/2)*i2+0];
870 if(verbose>5)
_printf_(
" writing Cracked vertices\n");
882 if(verbose>5)
_printf_(
" writing Cracked vertices\n");
894 xDelete<int>(connectivitysize_1);
895 xDelete<int>(head_1);
896 xDelete<int>(next_1);
897 xDelete<int>(connectivitysize_2);
898 xDelete<int>(head_2);
899 xDelete<int>(next_2);
◆ ReadMetric()
void bamg::Mesh::ReadMetric |
( |
const BamgOpts * |
bamgopts | ) |
|
Definition at line 904 of file Mesh.cpp.
916 if(verbose>3)
_printf_(
" processing metric\n");
919 double coef = bamgopts->
coeff;
934 a=bamgopts->
metric[i*3+0];
935 b=bamgopts->
metric[i*3+1];
936 c=bamgopts->
metric[i*3+2];
938 EigenMetric Vp(M/coef);
◆ WriteMetric()
void bamg::Mesh::WriteMetric |
( |
BamgOpts * |
bamgopts | ) |
|
◆ WriteIndex()
void bamg::Mesh::WriteIndex |
( |
int ** |
pindex, |
|
|
int * |
pnels |
|
) |
| |
Definition at line 960 of file Mesh.cpp.
977 index=xNew<int>(3*k);
980 if (t.det>0 && !(t.Hidden(0)||t.Hidden(1) || t.Hidden(2) )){
983 index[num*3+0]=
GetId(t[0])+1;
984 index[num*3+1]=
GetId(t[1])+1;
985 index[num*3+2]=
GetId(t[2])+1;
◆ AddMetric()
void bamg::Mesh::AddMetric |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 999 of file Mesh.cpp.
1009 if (Hessiantype==0){
1012 else if (Hessiantype==1){
1016 _error_(
"Hessiantype " << Hessiantype <<
" not supported yet (1->use Green formula, 0-> double L2 projection)");
◆ BuildMetric0()
void bamg::Mesh::BuildMetric0 |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 1638 of file Mesh.cpp.
1662 double* ss=(
double*)s;
1664 double* detT =
new double[
nbt];
1665 double* sumareas =
new double[
nbv];
1667 double* beta =
new double[
nbt*3];
1668 double* dx_elem =
new double[
nbt];
1669 double* dy_elem =
new double[
nbt];
1670 double* dx_vertex =
new double[
nbv];
1671 double* dy_vertex =
new double[
nbv];
1672 double* dxdx_elem =
new double[
nbt];
1673 double* dxdy_elem =
new double[
nbt];
1674 double* dydy_elem =
new double[
nbt];
1675 double* dxdx_vertex=
new double[
nbv];
1676 double* dxdy_vertex=
new double[
nbv];
1677 double* dydy_vertex=
new double[
nbv];
1681 _printf_(
" Construction of Metric: number of field: " << nbsol <<
" (nbt=" <<
nbt <<
", nbv=" <<
nbv <<
")\n");
1686 head_s=xNew<int>(
nbv);
1688 next_p=xNew<int>(3*
nbt);
1720 alpha[i*3+0]=(B.y-C.y)/dett;
1721 alpha[i*3+1]=(C.y-A.y)/dett;
1722 alpha[i*3+2]=(A.y-B.y)/dett;
1723 beta[ i*3+0]=(C.x-B.x)/dett;
1724 beta[ i*3+1]=(A.x-C.x)/dett;
1725 beta[ i*3+2]=(B.x-A.x)/dett;
1730 next_p[p]=head_s[k];
1741 for (
int nusol=0;nusol<nbsol;nusol++) {
1742 double smin=ss[nusol],smax=ss[nusol];
1745 for ( iv=0,k=0; iv<
nbv; iv++){
1746 smin=
Min(smin,ss[iv*nbsol+nusol]);
1747 smax=
Max(smax,ss[iv*nbsol+nusol]);
1749 double sdelta=smax-smin;
1753 if(verbose>2)
_printf_(
" Solution " << nusol <<
", Min = " << smin <<
", Max = " << smax <<
", Delta = " << sdelta <<
"\n");
1756 if (sdelta < 1.0e-10*
Max(absmax,1e-20)){
1757 _printf_(
" Solution " << nusol <<
" is constant, skipping...\n");
1762 for ( iv=0,k=0; iv<
nbv; iv++) dxdx_vertex[iv]=dxdy_vertex[iv]=dydy_vertex[iv]=dx_vertex[iv]=dy_vertex[iv]=0;
1765 for (i=0;i<
nbt;i++){
1773 sA = ss[iA*nbsol+nusol];
1774 sB = ss[iB*nbsol+nusol];
1775 sC = ss[iC*nbsol+nusol];
1779 dy_elem[i]=sA*beta[ 3*i+0]+sB*beta[ 3*i+1]+sC*beta[ 3*i+2];
1785 for(p=head_s[i];p!=-1;p=next_p[p]){
1788 dx_vertex[i]+=dx_elem[k]*detT[k]/sumareas[i];
1789 dy_vertex[i]+=dy_elem[k]*detT[k]/sumareas[i];
1794 for (i=0;i<
nbt;i++){
1802 dxdx_elem[i]=dx_vertex[iA]*
alpha[3*i+0]+dx_vertex[iB]*
alpha[3*i+1]+dx_vertex[iC]*
alpha[3*i+2];
1803 dxdy_elem[i]=dy_vertex[iA]*
alpha[3*i+0]+dy_vertex[iB]*
alpha[3*i+1]+dy_vertex[iC]*
alpha[3*i+2];
1804 dydy_elem[i]=dy_vertex[iA]*beta[3*i+0]+dy_vertex[iB]*beta[3*i+1]+dy_vertex[iC]*beta[3*i+2];
1810 for(p=head_s[i];p!=-1;p=next_p[p]){
1813 dxdx_vertex[i]+=dxdx_elem[k]*detT[k]/sumareas[i];
1814 dxdy_vertex[i]+=dxdy_elem[k]*detT[k]/sumareas[i];
1815 dydy_vertex[i]+=dydy_elem[k]*detT[k]/sumareas[i];
1820 for ( iv=0;iv<
nbv;iv++){
1821 vertices[iv].
MetricFromHessian(dxdx_vertex[iv],dxdy_vertex[iv],dydy_vertex[iv],smin,smax,ss[iv*nbsol+nusol],bamgopts->
err[nusol],bamgopts);
1827 xDelete<int>(head_s);
1828 xDelete<int>(next_p);
1835 delete [] dx_vertex;
1836 delete [] dy_vertex;
1837 delete [] dxdx_elem;
1838 delete [] dxdy_elem;
1839 delete [] dydy_elem;
1840 delete [] dxdx_vertex;
1841 delete [] dxdy_vertex;
1842 delete [] dydy_vertex;
◆ BuildMetric1()
void bamg::Mesh::BuildMetric1 |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 1845 of file Mesh.cpp.
1869 long i,k,iA,iB,iC,iv;
1871 double* ss=(
double*)s;
1873 double* detT =
new double[
nbt];
1874 double* Mmass=
new double[
nbv];
1875 double* Mmassxx=
new double[
nbv];
1876 double* dxdx=
new double[
nbv];
1877 double* dxdy=
new double[
nbv];
1878 double* dydy=
new double[
nbv];
1879 double* workT=
new double[
nbt];
1880 double* workV=
new double[
nbv];
1881 int* OnBoundary =
new int[
nbv];
1885 _printf_(
" Construction of Metric: number of field: " << nbsol <<
" (nbt=" <<
nbt <<
", nbv=" <<
nbv <<
")\n");
1889 for (iv=0;iv<
nbv;iv++){
1896 for (i=0;i<
nbt;i++){
1921 for(
int j=0;j<3;j++){
1923 Triangle *ta=t.Adj(j);
1925 if ( !ta || !ta->link){
1942 Mmassxx[iA] += dett;
1943 Mmassxx[iB] += dett;
1944 Mmassxx[iC] += dett;
1952 for (
int nusol=0;nusol<nbsol;nusol++) {
1954 double smin=ss[nusol],smax=ss[nusol];
1955 double h1=1.e30,h2=1e-30,rx=0;
1956 double hn1=1.e30,hn2=1e-30,rnx =1.e-30;
1959 for ( iv=0,k=0; iv<
nbv; iv++ ){
1960 dxdx[iv]=dxdy[iv]=dydy[iv]=0;
1961 smin=
Min(smin,ss[iv*nbsol+nusol]);
1962 smax=
Max(smax,ss[iv*nbsol+nusol]);
1964 double sdelta=smax-smin;
1968 if(verbose>2)
_printf_(
" Solution " << nusol <<
", Min = " << smin <<
", Max = " << smax <<
", Delta = " << sdelta <<
", number of fields = " << nbsol <<
"\n");
1971 if (sdelta < 1.0e-10*
Max(absmax,1e-20) ){
1972 if (verbose>2)
_printf_(
" Solution " << nusol <<
" is constant, skipping...\n");
1980 for ( iv=0,k=0; iv<
nbv; iv++) dxdx[iv]=dxdy[iv]=dydy[iv]=0;
1983 for (i=0;i<
nbt;i++){
2011 sA = ss[iA*nbsol+nusol];
2012 sB = ss[iB*nbsol+nusol];
2013 sC = ss[iC*nbsol+nusol];
2029 R2 Grads=(nAB*sC+nBC*sA+nCA*sB)/detT[i];
2034 if ( !tBC || !tBC->link ) nBC=O;
2035 if ( !tCA || !tCA->link ) nCA=O;
2036 if ( !tAB || !tAB->link ) nAB=O;
2042 dxdx[iA] += ( nCA.x + nAB.x ) *Grads.
x;
2043 dxdx[iB] += ( nAB.x + nBC.x ) *Grads.x;
2044 dxdx[iC] += ( nBC.x + nCA.x ) *Grads.x;
2047 dxdy[iA] += (( nCA.y + nAB.y ) *Grads.x + ( nCA.x + nAB.x ) *Grads.y) ;
2048 dxdy[iB] += (( nAB.y + nBC.y ) *Grads.x + ( nAB.x + nBC.x ) *Grads.y) ;
2049 dxdy[iC] += (( nBC.y + nCA.y ) *Grads.x + ( nBC.x + nCA.x ) *Grads.y) ;
2051 dydy[iA] += ( nCA.y + nAB.y ) *Grads.y;
2052 dydy[iB] += ( nAB.y + nBC.y ) *Grads.y;
2053 dydy[iC] += ( nBC.y + nCA.y ) *Grads.y;
2059 for ( iv=0,k=0 ; iv<
nbv; iv++){
2061 dxdx[iv] /= 2*Mmassxx[iv];
2063 dxdy[iv] /= 4*Mmassxx[iv];
2064 dydy[iv] /= 2*Mmassxx[iv];
2066 Metric M(dxdx[iv], dxdy[iv], dydy[iv]);
2080 for (
int xy = 0;xy<3;xy++) {
2082 else if (xy==1) dd=dxdy;
2083 else if (xy==2) dd=dydy;
2090 delete [] OnBoundary;
2094 for (
int ijacobi=0;ijacobi<
Max(NbJacobi,2);ijacobi++){
2103 cc =
Max((
double) ((Mmassxx[iA]>0)+(Mmassxx[iB]>0)+(Mmassxx[iC]>0)),1.);
2104 workT[i] = (dd[iA]+dd[iB]+dd[iC])/cc;
2106 for (iv=0;iv<
nbv;iv++) workV[iv]=0;
2108 for (i=0;i<
nbt;i++){
2114 double cc = workT[i]*detT[i];
2121 for (iv=0;iv<
nbv;iv++){
2122 if( ijacobi<NbJacobi || OnBoundary[iv]){
2123 dd[iv] = workV[iv]/(Mmass[iv]*6);
2130 for ( iv=0;iv<
nbv;iv++){
2144 delete [] OnBoundary;
◆ BuildGeometryFromMesh()
void bamg::Mesh::BuildGeometryFromMesh |
( |
BamgOpts * |
bamgopts = NULL | ) |
|
Definition at line 1197 of file Mesh.cpp.
1207 if(bamgopts) verbose=bamgopts->
verbose;
1210 if (verbose>1)
_printf_(
" construction of the geometry from the 2d mesh\n");
1213 if(
nbt<=0 ||
nbv <=0 )
_error_(
"nbt or nbv is negative (Mesh empty?)");
1218 SetOfEdges4* edge4=
new SetOfEdges4(
nbt*3,
nbv);
1219 long* st =
new long[
nbt*3];
1222 for (i=0;i<
nbt*3;i++) st[i]=-1;
1225 for (i=0;i<
nbe;i++){
1229 if (
nbe != edge4->nb()){
1231 _error_(
"Some Double edge in the mesh, the number is " <<
nbe <<
", nbe4=" << edge4->nb());
1237 for (i=0;i<
nbt;i++){
1245 if(st[k]==-1) st[k]=3*i+j;
1251 _error_(
"problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)");
1266 _printf_(
"Edge " << j <<
" of triangle " << i <<
"\n");
1267 _printf_(
"Edge " << (-st[k]+2)%3 <<
" of triangle " << (-st[k]+2)/3 <<
"\n");
1268 _printf_(
"Edge " <<
triangles[(-st[k]+2)/3].NuEdgeTriangleAdj((
int)((-st[k]+2)%3)) <<
" of triangle " <<
GetId(
triangles[(-st[k]+2)/3].TriangleAdj((
int)((-st[k]+2)%3))) <<
"\n");
1269 _error_(
"An edge belongs to more than 2 triangles");
1275 long nbedges = edge4->nb();
1276 delete edge4; edge4=NULL;
1281 _printf_(
" - number of vertices = " <<
nbv <<
"\n");
1282 _printf_(
" - number of triangles = " <<
nbt <<
"\n");
1283 _printf_(
" - number of given edges = " <<
nbe <<
"\n");
1284 _printf_(
" - number of all edges = " << nbedges <<
"\n");
1285 _printf_(
" - Euler number 1 - nb of holes = " <<
nbt-nbedges+
nbv <<
"\n");
1290 for (i=0;i<nbedges;i++){
1296 j = (int) ((-2-st[i])%3);
1310 Edge* edgessave=
edges;
1315 if(verbose>4)
_printf_(
" Construction of the edges " <<
nbe <<
"\n");
1317 for (i=0;i<nbedges;i++){
1323 j = (int) ((-2-st[i])%3);
1328 else if (st[i] >=0){
1330 j = (int) (st[i]%3);
1333 if (add>=0 && add <
nbe){
1350 _error_(
"problem in edge construction process: k!=nbe (should not happen)");
1353 if (edgessave)
delete [] edgessave;
1363 for (i=0;i<
nbe;i++){
1364 for (j=0;j<2;j++)
edges[i].v[j]->color++;
1368 for (i=0;i<
nbv;i++) {
1373 for (i=0;i<
nbe;i++){
1376 BamgVertex* v=
edges[i].
v[j];
1385 if(i0==-1) v->color=i*2+j;
1394 if (v!=
edges[i0 ].v[j0]){
1395 _error_(
"v!=edges[i0 ].v[j0]: this should not happen as the vertex belongs to this edge");
1415 long* colorT=
new long[
nbt];
1419 for (it=0;it<
nbt;it++) colorT[it]=-1;
1422 for (it=0;it<
nbt;it++){
1437 if( (j=st[level]++)<3 ){
1442 if ( ! t->Locked(j) && tt && (colorT[jt =
GetId(tt)] == -1) && ( tt->color==kolor)) {
1464 for (it=0;it<
nbt;it++){
1465 for (
int j=0;j<3;j++){
1488 long* colorV =
new long[
nbv];
1489 for (i=0;i<
nbv;i++) colorV[i]=-1;
1490 for (i=0;i<
nbe;i++){
1491 for ( j=0;j<2;j++) colorV[
GetId(
edges[i][j])]=0;
1495 for (i=0;i<
nbv;i++){
1496 if(!colorV[i]) colorV[i]=k++;
1506 if (verbose>3)
_printf_(
" number of vertices = " <<
Gh.
nbv <<
"\n number of edges = " <<
Gh.
nbe <<
"\n");
1513 for (i=0;i<
nbv;i++){
1514 if((j=colorV[i])>=0){
1526 for (i=0;i<
Gh.
nbv;i++) {
1537 long MaxICoord = 1073741823;
1541 _error_(
"Gh.coefIcoor<=0 in infered Geometry (this should not happen)");
1547 double * len =
new double[
Gh.
nbv];
1548 for(i=0;i<
Gh.
nbv;i++) len[i]=0;
1551 edge4=
new SetOfEdges4(
nbe,
nbv);
1552 double hmin = HUGE_VAL;
1554 for (i=0;i<
nbe;i++){
1558 long j0 = colorV[i0];
1559 long j1 = colorV[i1];
1570 if(required) kreq++;
1580 hmin =
Min(hmin,l12);
1587 hmin =
Min(hmin,l12);
1590 k = edge4->SortAndAdd(i0,i1);
1594 _error_(
"problem in Edge4 construction: k != i");
1599 for (i=0;i<
Gh.
nbv;i++){
1614 k = edge4->SortAndFind(i0,i1);
1617 _error_(
"This should not happen");
1630 for (i=0;i<
nbt;i++){
◆ RandomNumber()
long bamg::Mesh::RandomNumber |
( |
long |
max | ) |
|
◆ ReconstructExistingMesh()
void bamg::Mesh::ReconstructExistingMesh |
( |
BamgOpts * |
bamgopts | ) |
|
Definition at line 3302 of file Mesh.cpp.
3317 if(bamgopts) verbose=bamgopts->
verbose;
3323 if(verbose>2)
_printf_(
" Reconstruct mesh of " <<
nbv <<
" vertices\n");
3336 SetOfEdges4* edge4=
new SetOfEdges4(
nbt*3,
nbv);
3337 for (i=0;i<
nbe;i++){
3341 _error_(
"There are " << kk-
nbe <<
" double edges in the mesh");
3345 long* st =
new long[
nbt*3];
3346 for (i=0;i<
nbt*3;i++) st[i]=-1;
3347 for (i=0;i<
nbt;i++){
3348 for (
int j=0;j<3;j++){
3356 if(st[k]==-1) st[k]=3*i+j;
3381 _printf_(
" - number of vertices = " <<
nbv <<
" \n");
3382 _printf_(
" - number of triangles = " <<
nbt <<
" \n");
3383 _printf_(
" - number of given edges = " <<
nbe <<
" \n");
3384 _printf_(
" - number of all edges = " << edge4->nb() <<
"\n");
3385 _printf_(
" - Euler number 1 - nb of holes = " <<
nbt-edge4->nb()+
nbv <<
"\n");
3390 for (i=0;i<edge4->nb();i++){
3393 long i0=edge4->i(i);
3395 long i1=edge4->j(i);
3402 _printf_(
"Lost boundary edges " << i <<
" : " << edge4->i(i) <<
" " << edge4->j(i) <<
"\n");
3405 _printf_(
"Other lost boundary edges not shown...\n");
3411 _error_(k <<
" boundary edges (from the geometry) are not defined as mesh edges");
3416 for (i=0;i<
nbv;i++){
3428 long Nbtriafillhole=2*nbvb;
3429 Triangle* triafillhole=
new Triangle[Nbtriafillhole];
3438 _error_(
"ReconstructExistingMesh: All the vertices are aligned");
3470 quadtree =
new BamgQuadtree(
this,0);
3476 for (
int icount=2; icount<nbvb; icount++) {
3482 NbSwap += vi->Optim(1,1);
3486 AdjacentTriangle ta(0,0);
3487 long nbloss = 0,knbe=0;
3488 for ( i = 0; i <
nbe; i++){
3498 _error_(
"we lost " << nbloss <<
" existing edges other " << knbe);
3505 for (i=0;i<
nbt;i++){
3508 for (
int j=0;j<3;j++){
3510 Triangle &tta = *(Triangle*)ta;
3516 BamgVertex *v1= ta.EdgeVertex(1);
3520 tta.SetAdj2(ja,savetriangles + st[k] / 3,(
int) (st[k]%3));
3527 long NbTfillHoll =0;
3528 for (i=0;i<
nbt;i++){
3530 triangles[i]=Triangle((BamgVertex *) NULL,(BamgVertex *) NULL,(BamgVertex *) NULL);
3537 _assert_(savenbt+NbTfillHoll<=savemaxnbt);
3540 for (i=0;i<
nbt;i++){
3543 savetriangles[savenbt].
link=0;
3550 for (i=0;i<savenbt;i++) {
3551 Triangle & ti = savetriangles[i];
3552 for (
int j=0;j<3;j++){
3553 Triangle * ta = ti.TriangleAdj(j);
3554 int aa = ti.NuEdgeTriangleAdj(j);
3555 int lck = ti.Locked(j);
3557 else if ( ta >=
triangles && ta < tmax){
3558 ta= savetriangles + ta->color;
3559 ti.SetAdj2(j,ta,aa);
3560 if(lck) ti.SetLocked(j);
3573 _error_(
"number of triangles edges alone = " << k);
3584 for (i=0;i<
nbe;i++){
3586 if(
edges[i].GeomEdgeHook){
3587 for(
int j=0;j<2;j++){
3589 if (!
edges[i].adj[j]){
3592 if(!
edges[i][j].GeomEdgeHook->IsRequiredVertex()){
3593 _printf_(
"ReconstructExistingMesh error message: problem with the edge number " << i+1 <<
": [" <<
GetId(
edges[i][0])+1 <<
" " <<
GetId(
edges[i][1])+1 <<
"]\n");
3595 if (
edges[i][j].GeomEdgeHook->OnGeomVertex())
3597 else if (
edges[i][j].GeomEdgeHook->OnGeomEdge())
3600 _printf_(
"Its pointer is " <<
edges[i][j].GeomEdgeHook <<
"\n");
3602 _printf_(
"This edge is on geometry and has no adjacent edge (open curve) and one of the tip is not required\n");
3603 _error_(
"See above (might be cryptic...)");
◆ CreateSingleVertexToTriangleConnectivity()
void bamg::Mesh::CreateSingleVertexToTriangleConnectivity |
( |
| ) |
|
|
inline |
Definition at line 121 of file Mesh.h.
123 for (
int i=0;i<
nbt;i++)
triangles[i].SetSingleVertexToTriangleConnectivity();
◆ UnMarkUnSwapTriangle()
void bamg::Mesh::UnMarkUnSwapTriangle |
( |
| ) |
|
|
inline |
Definition at line 125 of file Mesh.h.
126 for (
int i=0;i<
nbt;i++)
◆ SetVertexFieldOn()
void bamg::Mesh::SetVertexFieldOn |
( |
| ) |
|
|
inline |
◆ SetVertexFieldOnBTh()
void bamg::Mesh::SetVertexFieldOnBTh |
( |
| ) |
|
|
inline |
◆ TriangulateFromGeom1()
void bamg::Mesh::TriangulateFromGeom1 |
( |
BamgOpts * |
bamgopts, |
|
|
int |
KeepVertices = 1 |
|
) |
| |
|
private |
Definition at line 4426 of file Mesh.cpp.
4433 if(bamgopts) verbose=bamgopts->
verbose;
4487 for (i=0;i<
Gh.
nbv;i++){
4488 if (
Gh[i].Required()) {
4495 else Gh[i].MeshVertexHook=0;
4499 if (vog.IsRequiredVertex()){
4501 BamgVertex *bv = vog;
4504 gv->MeshVertexHook->m = bv->m;
4514 for (
int i=0;i<
Gh.
nbcurves;i++) bcurve[i]=-1;
4517 for (
int iedge=0;iedge<
BTh.
nbe;iedge++){
4521 for(
int je=0;je<2;je++){
4524 if (ei[je].GeomEdgeHook->IsRequiredVertex()){
4538 bcurve[nc]=iedge*2+je;
4542 bcurve[nc]=iedge*2+je;
4550 _error_(
"problem generating number of curves (" <<
Gh.
nbcurves <<
" found in the geometry but " << bfind <<
" curve found in the mesh)");
4558 long nbex=0,NbVerticesOnGeomEdgex=0;
4559 for (
int step=0; step <2;step++){
4561 long NbOfNewPoints=0;
4568 for (
int icurve=0;icurve<
Gh.
nbcurves;icurve++){
4571 iedge=bcurve[icurve]/2;
4572 int jedge=bcurve[icurve]%2;
4579 long NbCreatePointOnCurve=0;
4582 for(
int phase=0;phase<=step;phase++){
4588 int icurveequi=
Gh.
GetId(curve);
4591 if(phase==0 && icurveequi!=icurve)
continue;
4595 int iedgeequi=bcurve[icurveequi]/2;
4596 int jedgeequi=bcurve[icurveequi]%2;
4598 int k0equi=jedgeequi,k1equi;
4599 Edge * peequi=
BTh.
edges+iedgeequi;
4605 GeomVertex * GA0 = *(*peequi)[k0equi].GeomEdgeHook;
4607 A0 = GA0->MeshVertexHook;
4610 Edge* PreviousNewEdge = 0;
4614 if(ongequi->Required()){
4615 GeomVertex *GA1 = *(*peequi)[1-k0equi].GeomEdgeHook;
4616 A1 = GA1->MeshVertexHook;
4621 Edge &eeequi=*peequi;
4625 ee.GeomEdgeHook->SetMark();
4626 BamgVertex & v0=ee[0], & v1=ee[1];
4634 while ((i!=NbCreatePointOnCurve) && sNew<=L) {
4647 double se= (sNew-L0)/LAB;
4648 if (se<0 || se>=1.000000001){
4649 _error_(
"Problem creating point on a boundary: se=" << se <<
" should be in [0 1]");
4653 _error_(
"Problem creating point on a boundary: se=" << se <<
" should be in [0 1]");
4655 se = k1 ? se : 1. - se;
4656 se = k1==k1equi ? se : 1. - se;
4660 e->GeomEdgeHook = ongequi;
4663 e->ReferenceNumber = eeequi.ReferenceNumber;
4664 e->adj[0]=PreviousNewEdge;
4666 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
4670 if (++i== NbCreatePointOnCurve)
break;
4675 _assert_(ee.GeomEdgeHook->CurveNumber==ei.GeomEdgeHook->CurveNumber);
4676 if (ee[k1].GeomEdgeHook->IsRequiredVertex()) {
4677 _assert_(eeequi[k1equi].GeomEdgeHook->IsRequiredVertex());
4678 GeomVertex * GA1 = *eeequi[k1equi].GeomEdgeHook;
4679 A1=GA1->MeshVertexHook;
4687 k0 = pe->Intersection(ee);
4688 peequi= eeequi.adj[k1equi];
4689 k0equi=peequi->Intersection(eeequi);
4695 e->GeomEdgeHook = ongequi;
4698 e->ReferenceNumber = peequi->ReferenceNumber;
4699 e->adj[0]=PreviousNewEdge;
4701 if (PreviousNewEdge) PreviousNewEdge->adj[1]=e;
4702 PreviousNewEdge = e;
4708 long NbSegOnCurve =
Max((
long)(L+0.5),(
long) 1);
4709 Lstep = L/NbSegOnCurve;
4710 NbCreatePointOnCurve = NbSegOnCurve-1;
4711 NbOfNewEdge += NbSegOnCurve;
4712 NbOfNewPoints += NbCreatePointOnCurve;
4720 _error_(
"too many vertices on geometry: " <<
nbv+NbOfNewPoints <<
" >= " <<
maxnbv);
4722 edges =
new Edge[NbOfNewEdge];
4728 NbVerticesOnGeomEdgex = NbOfNewPoints;
4738 if (verbose>4)
_printf_(
" -- current number of vertices = " <<
nbv <<
"\n");
4739 if (verbose>3)
_printf_(
" Creating initial Constrained Delaunay Triangulation...\n");
4740 if (verbose>3)
_printf_(
" Inserting boundary points\n");
4744 if (verbose>3)
_printf_(
" Forcing boundaries\n");
4748 if (verbose>3)
_printf_(
" Extracting subdomains\n");
4751 if (verbose>3)
_printf_(
" Inserting internal points\n");
4753 if (verbose>4)
_printf_(
" -- current number of vertices = " <<
nbv <<
"\n");
◆ TriangulateFromGeom0()
void bamg::Mesh::TriangulateFromGeom0 |
( |
BamgOpts * |
bamgopts | ) |
|
|
private |
Definition at line 4112 of file Mesh.cpp.
4120 int NbNewPoints,NbEdgeCurve;
4121 double lcurve,lstep,s;
4122 const int MaxSubEdge = 10;
4133 if(bamgopts) verbose=bamgopts->
verbose;
4136 bool background=(&
BTh !=
this);
4141 for (i=0;i<
Gh.
nbv;i++){
4149 for (i=0;i<
Gh.
nbv;i++){
4151 if (
Gh[i].Required()) {
4179 for (
int step=0;step<2;step++){
4189 for (i=0;i<
Gh.
nbe;i++){
4195 for(
int j=0;j<2;j++) {
4198 if (!ei.Mark() && ei[j].Required()){
4201 Edge* PreviousNewEdge=NULL;
4229 for (
int kstep=0;kstep<=step;kstep++){
4235 if (nbvend>=
maxnbv)
_error_(
"maximum number of vertices too low! Check the domain outline or increase maxnbv");
4243 va = a->MeshVertexHook;
4251 Metric MA = background ?
BTh.
MetricAt(a->r) :a->m ;
4252 Metric MB = background ?
BTh.
MetricAt(b->r) :b->m ;
4253 double ledge = (MA.Length(AB.x,AB.y) + MB.Length(AB.x,AB.y))/2;
4263 double lSubEdge[MaxSubEdge];
4268 lSubEdge[0] = ledge;
4273 NbSubEdge =
Min( MaxSubEdge, (
int) (ledge +0.5));
4283 double x =0, xstep= 1./NbSubEdge;
4284 for (
int kk=0; kk < NbSubEdge; kk++,A=B,MAs=MBs ) {
4286 B = e->F(k ? x : 1-x);
4287 MBs= background ?
BTh.
MetricAt(B) : Metric(1-x,MA,x,MB);
4289 lSubEdge[kk]=(ledge+=(MAs.Length(AB.x,AB.y)+MBs.Length(AB.x,AB.y))/2);
4293 double lcurveb = lcurve+ledge;
4296 while(s>=lcurve && s<=lcurveb &&
nbv<nbvend){
4307 double ss = s-lcurve;
4310 int kk0=-1,kk1=NbSubEdge-1,kkk;
4311 double ll0=0,ll1=ledge,llk;
4313 if (ss < (llk=lSubEdge[kkk=(kk0+kk1)/2]))
4325 double sbb = (ss-ll0)/(ll1-ll0);
4327 double bb = (kk1+sbb)/NbSubEdge;
4332 vb->
m = Metric(aa,a->m,bb,b->m);
4333 vb->ReferenceNumber = e->ReferenceNumber;
4334 double abcisse = k ? bb : aa;
4335 vb->r = e->F(abcisse);
4345 if(PreviousNewEdge) PreviousNewEdge->
adj[1]=&
edges[
nbe];
4357 if (b->Required())
break;
4359 k = e->AdjVertexIndex[kprev];
4363 vb = b->MeshVertexHook;
4366 NbEdgeCurve =
Max((
long) (lcurve +0.5), (
long) 1);
4368 NbNewPoints = NbEdgeCurve-1;
4370 NbVerticesOnGeomEdge0 += NbNewPoints;
4373 nbvend=
nbv+NbNewPoints;
4374 lstep = lcurve / NbEdgeCurve;
4383 if(PreviousNewEdge) PreviousNewEdge->
adj[1] = &
edges[
nbe];
4386 else nbe += NbEdgeCurve;
4396 if(NbVerticesOnGeomEdge0)
VerticesOnGeomEdge =
new VertexOnGeom[NbVerticesOnGeomEdge0];
4408 if (verbose>4)
_printf_(
" -- current number of vertices = " <<
nbv <<
"\n");
4409 if (verbose>3)
_printf_(
" Creating initial Constrained Delaunay Triangulation...\n");
4410 if (verbose>3)
_printf_(
" Inserting boundary points\n");
4414 if (verbose>3)
_printf_(
" Forcing boundaries\n");
4418 if (verbose>3)
_printf_(
" Extracting subdomains\n");
4421 if (verbose>3)
_printf_(
" Inserting internal points\n");
4423 if (verbose>4)
_printf_(
" -- current number of vertices = " <<
nbv <<
"\n");
◆ Triangulate()
void bamg::Mesh::Triangulate |
( |
double * |
x, |
|
|
double * |
y, |
|
|
int |
nods, |
|
|
BamgOpts * |
bamgopts |
|
) |
| |
|
private |
Definition at line 4084 of file Mesh.cpp.
4091 if(bamgopts) verbose=bamgopts->
verbose;
4098 if(verbose)
_printf_(
"Reading vertices (" <<
nbv <<
")\n");
◆ Init()
void bamg::Mesh::Init |
( |
long |
maxnbv_in | ) |
|
|
private |
◆ ForceEdge()
Definition at line 4763 of file Mesh.cpp.
4771 taret=AdjacentTriangle(0,0);
4774 BamgVertex *v1, *v2 = tta.EdgeVertex(0),*vbegin =v2;
4777 long long det2 = v2 ?
det(*v2,a,b): -1 , det1;
4779 det2 =
det(*v2,a,b);
4782 v2 = tta.EdgeVertex(0);
4787 det2 =
det(*v2,a,b);
4793 v2 = tc.EdgeVertex(0);
4795 det2 = v2 ?
det(*v2,a,b): det2;
4797 if((det1 < 0) && (det2 >0)) {
4799 BamgVertex * va = &a, *vb = &b;
4804 long long detss = 0,l=0;
4806 if(l++ > 10000000) {
4807 _error_(
"Loop in forcing Egde, nb de swap=" << NbSwap <<
", nb of try swap (" << l <<
") too big");
4809 BamgVertex *aa = tc.EdgeVertex(0), *bb = tc.EdgeVertex(1);
4810 if (((aa == &a ) && (bb == &b)) ||((bb == &a ) && (aa == &b))){
4828 if ( vbegin == v2 )
return -1;
◆ SwapForForcingEdge()
Definition at line 4838 of file Mesh.cpp.
4846 if(tt1.Locked())
return 0;
4848 AdjacentTriangle tt2 =
Adj(tt1);
4849 Triangle *t1=tt1,*t2=tt2;
4850 short a1=tt1,a2=tt2;
4851 if ( a1<0 || a1>=3 ){
4859 long long dets2 =
det(*pva,*pvb,s2);
4860 long long det1=t1->det , det2=t2->det ;
4861 long long detT = det1+det2;
4862 if ((det1<=0 ) || (det2<=0)){
4863 _error_(
"(det1<=0 ) || (det2<=0)");
4865 if ( (detsa>=0) || (detsb<=0) ){
4866 _error_(
"(detsa>=0) || (detsb<=0)");
4869 long long ndet2 = detT - ndet1;
4872 if ((ndet1 >0) && (ndet2 >0))
4874 if ((dets1 <=0 && dets2 <=0) || (dets2 >=0 && detsb >=0))
4880 if (ToSwap) NbSwap++,
4886 dets1 = ToSwap ? dets1 : detsa ;
4889 else if (dets2 > 0){
4890 dets1 = ToSwap ? dets1 : detsb ;
4893 if(!ToSwap) tt1 =
Next(tt2);
4908 dets1 = (ToSwap ? dets1 : detsa) ;
4912 dets1 = (ToSwap ? dets1 : detsb) ;
4914 if(!ToSwap) tt1 =
Next(tt2);
◆ Gh
◆ BTh
◆ vertices
◆ triangles
◆ edges
◆ quadtree
◆ orderedvertices
◆ subdomains
◆ NbRef
◆ maxnbv
◆ maxnbt
◆ nbv
◆ nbt
◆ nbe
◆ nbsubdomains
long bamg::Mesh::nbsubdomains |
◆ nbtout
◆ pmin
◆ pmax
◆ coefIcoor
double bamg::Mesh::coefIcoor |
◆ lIntTria
◆ randomseed
long bamg::Mesh::randomseed |
◆ NbVerticesOnGeomVertex
long bamg::Mesh::NbVerticesOnGeomVertex |
◆ VerticesOnGeomVertex
◆ NbVerticesOnGeomEdge
long bamg::Mesh::NbVerticesOnGeomEdge |
◆ VerticesOnGeomEdge
◆ NbVertexOnBThVertex
long bamg::Mesh::NbVertexOnBThVertex |
◆ VertexOnBThVertex
◆ NbVertexOnBThEdge
long bamg::Mesh::NbVertexOnBThEdge |
◆ VertexOnBThEdge
◆ NbCrackedVertices
long bamg::Mesh::NbCrackedVertices |
◆ CrackedVertices
long* bamg::Mesh::CrackedVertices |
◆ NbCrackedEdges
long bamg::Mesh::NbCrackedEdges |
◆ CrackedEdges
The documentation for this class was generated from the following files:
long GetId(const Triangle &t) const
Edge ** MakeGeomEdgeToEdge()
void MetricFromHessian(const double Hxx, const double Hyx, const double Hyy, const double smin, const double smax, const double s, const double err, BamgOpts *bamgopts)
long NbVerticesOnGeomVertex
BamgVertex * MeshVertexHook
int VerticesOnGeomVertexSize[2]
void SetSingleVertexToTriangleConnectivity()
static const short PreviousEdge[3]
void Triangulate(double *x, double *y, int nods, BamgOpts *bamgopts)
#define _printf_(StreamArgs)
int GetReferenceNumber() const
int NodalConnectivitySize[2]
double * NodalElementConnectivity
double abscisseInterpole(const Metric &Ma, const Metric &Mb, R2 AB, double s, int optim)
int VerticesOnGeomEdgeSize[2]
long RandomNumber(long max)
double * NodalConnectivity
AdjacentTriangle Adj() const
int ForceEdge(BamgVertex &a, BamgVertex &b, AdjacentTriangle &taret)
void ReadMesh(int *index, double *x, double *y, int nods, int nels, BamgOpts *bamgopts)
AdjacentTriangle Adj(const AdjacentTriangle &a)
static const short PreviousVertex[3]
CrackedEdge * CrackedEdges
VertexOnEdge * VertexOnBThEdge
Metric MetricAt(const R2 &)
void AddVertex(BamgVertex &s, Triangle *t, long long *=0)
void ForceBoundary(BamgOpts *bamgopts)
VertexOnGeom * VerticesOnGeomVertex
static const short NextVertex[3]
ListofIntersectionTriangles lIntTria
double * ElementConnectivity
void Exchange(T &a, T &b)
BamgVertex * TooClose(BamgVertex *, double, int, int)
long BigPrimeNumber(long n)
T Max(const T &a, const T &b)
int EdgesOnGeomEdgeSize[2]
GeomSubDomain * subdomains
void SetIntCoor(const char *from=0)
void TriangulateFromGeom0(BamgOpts *bamgopts)
void BuildMetric1(BamgOpts *bamgopts)
int CrackedVerticesSize[2]
static const short OppositeEdge[3]
void BuildGeometryFromMesh(BamgOpts *bamgopts=NULL)
double * VerticesOnGeomEdge
GeomEdge * ProjectOnCurve(const Edge &, double, BamgVertex &, VertexOnGeom &) const
AdjacentTriangle Adj(int i) const
void ReconstructExistingMesh(BamgOpts *bamgopts)
RR Area2(const P2< R, RR > a, const P2< R, RR > b, const P2< R, RR > c)
void PostRead(bool checkcurve=false)
AdjacentTriangle Next(const AdjacentTriangle &ta)
double * SubDomainsFromGeom
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
void ReadGeometry(BamgGeom *bamggeom, BamgOpts *bamgopts)
void CreateSingleVertexToTriangleConnectivity()
long long det(const I2 &a, const I2 &b, const I2 &c)
RR Norme2(const P2< R, RR > x)
long GetId(const GeomVertex &t) const
void BuildMetric0(BamgOpts *bamgopts)
double * PreviousNumbering
long NewPoints(BamgVertex *, long &nbv, long maxnbv)
double LengthInterpole(const Metric &Ma, const Metric &Mb, R2 AB)
int SubDomainsFromGeomSize[2]
const int IsVertexOnVertex
void NewPoints(Mesh &, BamgOpts *bamgopts, int KeepVertices=1)
BamgVertex * NearestVertex(int i, int j)
void TriangleIntNumbering(long *renumbering)
static const short EdgesVertexTriangle[3][2]
double * VerticesOnGeomVertex
#define _error_(StreamArgs)
R2 I2ToR2(const I2 &P) const
AdjacentTriangle CloseBoundaryEdge(I2 A, Triangle *t, double &a, double &b)
static const short VerticesOfTriangularEdge[3][2]
I2 R2ToI2(const R2 &P) const
long NbVerticesOnGeomEdge
int SwapForForcingEdge(BamgVertex *&pva, BamgVertex *&pvb, AdjacentTriangle &tt1, long long &dets1, long long &detsa, long long &detsb, int &nbswap)
void FindSubDomain(BamgOpts *bamgopts, int OutSide=0)
long TriangleReferenceList(long *) const
int NodalElementConnectivitySize[2]
BamgVertex * EdgeVertex(const int &i) const
P2< R, RR > Orthogonal(const P2< R, RR > x)
void SplitEdge(Mesh &, const R2 &, const R2 &, int nbegin=0)
long InsertNewPoints(long nbvold, long &NbTSwap, BamgOpts *bamgopts)
Triangle * TriangleAdj(int i) const
Triangle * TriangleFindFromCoord(const I2 &, long long[3], Triangle *tstart=0)
void swap(Triangle *t1, short a1, Triangle *t2, short a2, BamgVertex *s1, BamgVertex *s2, long long det1, long long det2)
P2xP2< double, double > D2xD2
IssmDouble max(IssmDouble a, IssmDouble b)
void TriangulateFromGeom1(BamgOpts *bamgopts, int KeepVertices=1)
AdjacentTriangle Previous(const AdjacentTriangle &ta)
BamgVertex ** orderedvertices
T Min(const T &a, const T &b)
int ElementConnectivitySize[2]
static const short OppositeVertex[3]
void SetAdj2(short a, Triangle *t, short aat)
void Insert(BamgOpts *bamgopts)
VertexOnVertex * VertexOnBThVertex
static const short NextEdge[3]
VertexOnGeom * VerticesOnGeomEdge