![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 5448
- Timestamp:
- 08/20/10 11:22:07 (15 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/c/objects/Bamg/Mesh.cpp (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Bamg/Mesh.cpp
r5447 r5448 4875 4875 /*Generate mesh from geometry*/ 4876 4876 4877 Gh.NbRef++;// add a ref to GH 4878 4879 int i,j,k; 4880 int nbcurves=0,NbNewPoints,NbEdgeCurve; 4881 double lcurve,lstep,s; 4882 const int MaxSubEdge = 10; 4883 4884 R2 AB; 4885 GeometricalVertex *a,*b; 4886 BamgVertex *va,*vb; 4887 GeometricalEdge *e; 4877 /*Intermediaries*/ 4878 int i,j,k; 4879 int nbcurves = 0; 4880 int NbNewPoints,NbEdgeCurve; 4881 double lcurve,lstep,s; 4882 const int MaxSubEdge = 10; 4883 4884 R2 AB; 4885 GeometricalVertex *a, *b; 4886 BamgVertex *va, *vb; 4887 GeometricalEdge *e; 4888 4889 // add a ref to GH to make sure that it is not destroyed by mistake 4890 Gh.NbRef++; 4888 4891 4889 4892 /*Get options*/ 4890 4893 int verbose=bamgopts->verbose; 4891 4894 4892 //initialize this 4893 if (verbose>3) printf(" Generating Boundary vertices\n"); 4895 //initialize Mesh 4894 4896 Init(imaxnbv); 4895 4897 nbv=0; … … 4898 4900 4899 4901 //build background mesh flag (1 if background, else 0) 4900 int background=(&BTh != this); 4901 4902 //Compute number of vertices on geometrical vertex 4902 bool background=(&BTh != this); 4903 4904 /*Build VerticesOnGeomVertex*/ 4905 4906 //Compute the number of geometrical vertices that we are going to use to mesh 4903 4907 for (i=0;i<Gh.nbv;i++){ 4904 4908 if (Gh[i].Required()) NbVerticesOnGeomVertex++; 4905 4909 } 4906 4907 //initialize VerticesOnGeomVertex 4910 //allocate 4908 4911 VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 4909 if( NbVerticesOnGeomVertex >= maxnbv) { 4910 ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv); 4911 } 4912 4913 //Add all the geometrical vertices to the mesh 4912 if(NbVerticesOnGeomVertex >= maxnbv) ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv); 4913 //Build VerticesOnGeomVertex 4914 4914 nbv=0; 4915 4915 for (i=0;i<Gh.nbv;i++){ … … 4917 4917 if (Gh[i].Required()) {//Gh vertices Required 4918 4918 4919 //Add the vertex (provided that nbv<maxnbv) 4920 if (nbv<maxnbv){ 4921 vertices[nbv]=Gh[i]; 4922 } 4923 else{ 4924 ISSMERROR("Maximum number of vertices (maxnbv = %i) too small",maxnbv); 4925 } 4919 //Add the vertex 4920 ISSMASSERT(nbv<maxnbv); 4921 vertices[nbv]=Gh[i]; 4926 4922 4927 4923 //Add pointer from geometry (Gh) to vertex from mesh (Th) … … 4929 4925 4930 4926 //Build VerticesOnGeomVertex for current point 4931 VerticesOnGeomVertex[nbv]= 4927 VerticesOnGeomVertex[nbv]=VertexOnGeom(*Gh[i].to,Gh[i]); 4932 4928 4933 4929 //nbv increment … … 4936 4932 } 4937 4933 4938 / /check that edges has been allocated4939 if (edges){ 4940 ISSMERROR("edges is empty");4941 }4934 /*Build VerticesOnGeomEdge*/ 4935 4936 //check that edges is still empty (Init) 4937 ISSMASSERT(!edges); 4942 4938 4943 4939 /* Now we are going to create the first edges corresponding … … 4957 4953 4958 4954 //go through the edges of the geometry 4959 for (i=0;i<Gh.nbe;i++) 4955 for (i=0;i<Gh.nbe;i++){ 4960 4956 4961 4957 //ei = current Geometrical edge 4962 4958 GeometricalEdge &ei=Gh.edges[i]; 4963 4959 4960 //loop over the two vertices of the edge ei 4964 4961 for(int j=0;j<2;j++) { 4965 4962 4966 /*T he first time, no edge is marked but this might change during the loop*/4963 /*Take only required vertices (corner)*/ 4967 4964 if (!ei.Mark() && ei[j].Required()){ 4968 4965 4969 4966 long nbvend=0; 4970 4967 Edge* PreviousNewEdge=NULL; 4971 4972 4968 lstep = -1; 4973 4969 4974 /*If Edge is required */4975 if(ei.Required() ){4970 /*If Edge is required (do that only once for the 2 vertices)*/ 4971 if(ei.Required() && j==0){ 4976 4972 //do not create internal points if required (take it as is) 4977 if (j==0){ 4978 if(step==0) nbe++; 4979 else{ 4980 e=&ei; 4981 a=ei(0); 4982 b=ei(1); 4983 4984 //check that edges has been allocated 4985 if (!edges) ISSMERROR("edges has not been allocated..."); 4986 edges[nbe].v[0]=a->to; 4987 edges[nbe].v[1]=b->to;; 4988 edges[nbe].ReferenceNumber = e->ReferenceNumber; 4989 edges[nbe].GeometricalEdgeHook = e; 4990 edges[nbe].adj[0] = 0; 4991 edges[nbe].adj[1] = 0; 4992 nbe++; 4993 } 4973 if(step==0) nbe++; 4974 else{ 4975 e=&ei; 4976 a=ei(0); 4977 b=ei(1); 4978 4979 //check that edges has been allocated 4980 ISSMASSERT(edges); 4981 edges[nbe].v[0]=a->to; 4982 edges[nbe].v[1]=b->to;; 4983 edges[nbe].ReferenceNumber = e->ReferenceNumber; 4984 edges[nbe].GeometricalEdgeHook = e; 4985 edges[nbe].adj[0] = 0; 4986 edges[nbe].adj[1] = 0; 4987 nbe++; 4994 4988 } 4995 4989 } 4996 4990 4997 /*If Edge is not required: on a curve*/4991 /*If Edge is not required: we are on a curve*/ 4998 4992 else { 4999 4993 for (int kstep=0;kstep<=step;kstep++){ 5000 // step=0, do nothing5001 // step=1, compute the length of the curve5002 // step=2 create the points and edge4994 //kstep=0, do nothing 4995 //kstep=1, compute the length of the curve 4996 //kstep=2 create the points and edge 5003 4997 PreviousNewEdge=0; 5004 4998 NbNewPoints=0; … … 5008 5002 s = lstep; 5009 5003 5010 // i = edge number, j=[0;1] vertex number in edge 5011 5012 k=j; // k = vertex number in edge (0 or 1) 5004 /*reminder: i = edge number, j=[0;1] vertex index in edge*/ 5005 k=j; // k = vertex index in edge (0 or 1) 5013 5006 e=&ei; // e = reference of current edge 5014 5007 a=ei(k); // a = pointer toward the kth vertex of the current edge 5015 va = a->to; // va = pointer toward vertex associated5008 va = a->to; // va = pointer toward mesh vertex associated 5016 5009 e->SetMark(); // Mark edge 5017 5010 5018 //if SameGeo We have go to the background geometry 5019 //to find the discretisation of the curve 5011 /*If we have a Background mesh, we can use it to discretize the curve*/ 5020 5012 for(;;){ 5021 k = 1-k; 5022 b = (*e)(k); // b = pointer toward the other vertex of the current edge5013 k = 1-k; // other vertx index of the curve 5014 b = (*e)(k); // b = pointer toward the other vertex of the current edge 5023 5015 AB= b->r - a->r; // AB = vector of the current edge 5024 5016 Metric MA = background ? BTh.MetricAt(a->r) :a->m ; //Get metric associated to A 5025 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 5026 double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 5027 5028 /* We are now creating the edges of the mesh from the 5029 * geometrical edge selected above. 5030 * The edge will be divided according to the metric 5031 * previously computed and cannot be divided more 5032 * than 10 times (MaxSubEdge). */ 5017 Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B 5018 double ledge = (MA(AB) + MB(AB))/2; //Get edge length in metric 5019 5020 /* We are now creating the mesh edges from the geometrical edge selected above. 5021 * The edge will be divided according to the metric previously computed and cannot 5022 * be divided more than 10 times (MaxSubEdge). */ 5033 5023 5034 5024 //By default, there is only one subedge that is the geometrical edge itself … … 5039 5029 5040 5030 //Build Subedges according to the edge length 5041 //if ledge < 1.5 (between one and 2), take the edge as is 5042 if (ledge < 1.5) lSubEdge[0] = ledge; 5031 if (ledge < 1.5){ 5032 //if ledge < 1.5 (between one and 2), take the edge as is 5033 lSubEdge[0] = ledge; 5034 } 5043 5035 //else, divide the edge 5044 5036 else { … … 5054 5046 x += xstep; 5055 5047 B = e->F(k ? x : 1-x); 5056 MBs= background ? BTh.MetricAt(B) : Metric(1-x, MA, x ,MB);5048 MBs= background ? BTh.MetricAt(B) : Metric(1-x, MA, x ,MB); 5057 5049 AB = A-B; 5058 5050 lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2); … … 5138 5130 } // for (i=0;i<nbe;i++) 5139 5131 if(!step) { 5140 if (edges){ 5141 ISSMERROR("edges"); 5142 } 5143 if (VerticesOnGeomEdge){ 5144 ISSMERROR("VerticesOnGeomEdge"); 5145 } 5132 ISSMASSERT(!edges); 5133 ISSMASSERT(!VerticesOnGeomEdge); 5134 5146 5135 edges = new Edge[nbex=nbe]; 5147 if(NbVerticesOnGeomEdge0) 5148 VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0]; 5149 if (!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0){ 5150 ISSMERROR("!VerticesOnGeomEdge && NbVerticesOnGeomEdge0!=0"); 5151 } 5136 if(NbVerticesOnGeomEdge0) VerticesOnGeomEdge = new VertexOnGeom[NbVerticesOnGeomEdge0]; 5137 5152 5138 // do the vertex on a geometrical vertex 5139 ISSMASSERT(VerticesOnGeomEdge || NbVerticesOnGeomEdge0==0); 5153 5140 NbVerticesOnGeomEdge0 = NbVerticesOnGeomEdge; 5154 5141 } 5155 else if (NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0){5156 ISSM ERROR("NbVerticesOnGeomEdge != NbVerticesOnGeomEdge0");5142 else{ 5143 ISSMASSERT(NbVerticesOnGeomEdge==NbVerticesOnGeomEdge0); 5157 5144 } 5158 5145 }
Note:
See TracChangeset
for help on using the changeset viewer.
