Changeset 5531


Ignore:
Timestamp:
08/24/10 08:40:17 (15 years ago)
Author:
Mathieu Morlighem
Message:

As usual (and fixed a bug introduced by me...)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Bamg/Mesh.cpp

    r5489 r5531  
    99
    1010        static const  Direction NoDirOfSearch=Direction();
    11         long NbUnSwap =0;
    1211
    1312        /*Constructors/Destructors*/
     
    48954894                int verbose=bamgopts->verbose;
    48964895
    4897                 //initialize Mesh
    4898                 nbv=0;
    4899                 NbVerticesOnGeomVertex=0;
    4900                 NbVerticesOnGeomEdge=0;
    4901 
    49024896                //build background mesh flag (1 if background, else 0)
    49034897                bool background=(&BTh != this);
     
    49124906                VerticesOnGeomVertex = new VertexOnGeom[NbVerticesOnGeomVertex]; 
    49134907                if(NbVerticesOnGeomVertex >= maxnbv) ISSMERROR("too many vertices on geometry: %i >= %i",NbVerticesOnGeomVertex,maxnbv);
     4908                ISSMASSERT(nbv==0);
    49144909                //Build VerticesOnGeomVertex
    4915                 nbv=0;
    49164910                for (i=0;i<Gh.nbv;i++){
    49174911                        /* Add vertex only if required*/
     
    49414935                 * to the one present in the geometry provided.
    49424936                 * We proceed in 2 steps
    4943                  *  -step 1: we count all the edges
    4944                  *           we allocate the number of edges at the end of step 1
    4945                  *  -step 2: the edges are created */
     4937                 *  -step 0: we count all the edges
     4938                 *           we allocate the number of edges at the end of step 0
     4939                 *  -step 1: the edges are created */
    49464940                for (int step=0;step<2;step++){
    49474941
     
    49624956                                for(int j=0;j<2;j++) {
    49634957
    4964                                         /*Take only required vertices (corner)*/
     4958                                        /*Take only required vertices (corner->beginning of a new curve)*/
    49654959                                        if (!ei.Mark() && ei[j].Required()){
    49664960
     
    49704964
    49714965                                                /*If Edge is required (do that only once for the 2 vertices)*/
    4972                                                 if(ei.Required() && j==0){
    4973                                                         //do not create internal points if required (take it as is)
    4974                                                         if(step==0) nbe++;
    4975                                                         else{
    4976                                                                 e=&ei;
    4977                                                                 a=ei(0);
    4978                                                                 b=ei(1);
    4979 
    4980                                                                 //check that edges has been allocated
    4981                                                                 ISSMASSERT(edges);
    4982                                                                 edges[nbe].v[0]=a->MeshVertexHook;
    4983                                                                 edges[nbe].v[1]=b->MeshVertexHook;;
    4984                                                                 edges[nbe].ReferenceNumber = e->ReferenceNumber;
    4985                                                                 edges[nbe].GeometricalEdgeHook = e;
    4986                                                                 edges[nbe].adj[0] = 0;
    4987                                                                 edges[nbe].adj[1] = 0;
    4988                                                                 nbe++;
     4966                                                if(ei.Required()){
     4967                                                        if (j==0){
     4968                                                                //do not create internal points if required (take it as is)
     4969                                                                if(step==0) nbe++;
     4970                                                                else{
     4971                                                                        e=&ei;
     4972                                                                        a=ei(0);
     4973                                                                        b=ei(1);
     4974
     4975                                                                        //check that edges has been allocated
     4976                                                                        ISSMASSERT(edges);
     4977                                                                        edges[nbe].v[0]=a->MeshVertexHook;
     4978                                                                        edges[nbe].v[1]=b->MeshVertexHook;;
     4979                                                                        edges[nbe].ReferenceNumber = e->ReferenceNumber;
     4980                                                                        edges[nbe].GeometricalEdgeHook = e;
     4981                                                                        edges[nbe].adj[0] = 0;
     4982                                                                        edges[nbe].adj[1] = 0;
     4983                                                                        nbe++;
     4984                                                                }
    49894985                                                        }
    49904986                                                }
     
    49934989                                                else {
    49944990                                                        for (int kstep=0;kstep<=step;kstep++){
    4995                                                                 //kstep=0, do nothing
    4996                                                                 //kstep=1, compute the length of the curve
    4997                                                                 //kstep=2  create the points and edge
     4991                                                                //kstep=0, compute number of edges (discretize curve)
     4992                                                                //kstep=1  create the points and edge
    49984993                                                                PreviousNewEdge=0;
    49994994                                                                NbNewPoints=0;
     
    50014996                                                                if (nbvend>=maxnbv) ISSMERROR("maximum number of vertices too low! Check the domain outline or increase maxnbv");
    50024997                                                                lcurve =0;
    5003                                                                 s = lstep;
     4998                                                                s = lstep; //-1 initially, then length of each sub edge
    50044999
    50055000                                                                /*reminder: i = edge number, j=[0;1] vertex index in edge*/
     
    50105005                                                                e->SetMark();   // Mark edge
    50115006
    5012                                                                 /*If we have a Background mesh, we can use it to discretize the curve*/
     5007                                                                /*Loop until we reach the end of the curve*/
    50135008                                                                for(;;){
    50145009                                                                        k = 1-k;            // other vertx index of the curve
     
    50165011                                                                        AB= b->r - a->r;   // AB = vector of the current edge
    50175012                                                                        Metric MA = background ? BTh.MetricAt(a->r) :a->m ;  //Get metric associated to A
    5018                                                                         Metric MB = background ? BTh.MetricAt(b->r) :b->m ; //Get metric associated to B
    5019                                                                         double ledge = (MA(AB) + MB(AB))/2;                 //Get edge length in metric
     5013                                                                        Metric MB = background ? BTh.MetricAt(b->r) :b->m ;  //Get metric associated to B
     5014                                                                        double ledge = (MA(AB) + MB(AB))/2;                  //Get edge length in metric
    50205015
    50215016                                                                        /* We are now creating the mesh edges from the geometrical edge selected above.
     
    50365031                                                                        //else, divide the edge
    50375032                                                                        else {
    5038                                                                                 //compute number of subedges (division of the edge)
     5033                                                                                //compute number of subedges (division of the edge), Maximum is 10
    50395034                                                                                NbSubEdge = Min( MaxSubEdge, (int) (ledge +0.5));
    5040                                                                                 //A and B are the position of points on the edge
     5035                                                                                /*Now, we are going to divide the edge according to the metric.
     5036                                                                                 * Get segment by sement along the edge.
     5037                                                                                 * Build lSubEdge, which holds the distance between the first vertex
     5038                                                                                 * of the edge and the next point on the edge according to the
     5039                                                                                 * discretization (each SubEdge is AB)*/
    50415040                                                                                R2 A,B;
    50425041                                                                                A=a->r;
     
    50495048                                                                                        MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
    50505049                                                                                        AB = A-B;
    5051                                                                                         lSubEdge[kk]= (ledge+=(MAs(AB)+MBs(AB))/2);
     5050                                                                                        lSubEdge[kk]=(ledge+=(MAs(AB)+MBs(AB))/2);
    50525051                                                                                }
    50535052                                                                        }
     
    50565055
    50575056                                                                        /*Now, create corresponding points*/
    5058                                                                         while (lcurve<=s && s<=lcurveb && nbv<nbvend){
     5057                                                                        while(s>=lcurve && s<=lcurveb && nbv<nbvend){
    50595058
    50605059                                                                                double ss = s-lcurve;
     
    50955094                                                                                va = vb;
    50965095                                                                        }
     5096                                                                       
     5097                                                                        /*We just added one edge to the curve: Go to the next one*/
    50975098                                                                        lcurve = lcurveb;
    50985099                                                                        e->SetMark();
    50995100                                                                        a=b;
    5100                                                                         if (b->Required() ) break;
     5101
     5102                                                                        /*If b is required, we are on a new curve->break*/
     5103                                                                        if (b->Required()) break;
    51015104                                                                        int kprev=k;
    51025105                                                                        k = e->AdjVertexIndex[kprev];// next vertices
     
    51055108                                                                }// for(;;)
    51065109                                                                vb = b->MeshVertexHook;
     5110
     5111                                                                /*Number of edges in the last disretized curve*/
    51075112                                                                NbEdgeCurve = Max((long) (lcurve +0.5), (long) 1);
     5113                                                                /*Number of internal vertices in the last disretized curve*/
    51085114                                                                NbNewPoints = NbEdgeCurve-1;
    51095115                                                                if(!kstep){
     
    51125118                                                                }
    51135119                                                                nbvend=nbv+NbNewPoints;
    5114                                                                 lstep = lcurve / NbEdgeCurve;
     5120                                                                lstep = lcurve / NbEdgeCurve; //approximately one
    51155121                                                        }// end of curve --
    51165122                                                        if (edges) { // last edges of the curves
Note: See TracChangeset for help on using the changeset viewer.