Changeset 15611


Ignore:
Timestamp:
07/25/13 09:32:50 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: make the distinction between faces and edges for 3d Quadratic elements (quadratic requires edges, DG requires faces, the format is not the same)

Location:
issm/trunk-jpl/src/c
Files:
2 added
1 deleted
14 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r15464 r15611  
    244244                                        ./modules/ModelProcessorx/NodesPartitioning.cpp\
    245245                                        ./modules/ModelProcessorx/EdgesPartitioning.cpp\
     246                                        ./modules/ModelProcessorx/FacesPartitioning.cpp\
    246247                                        ./modules/ModelProcessorx/SortDataSets.cpp\
    247248                                        ./modules/ModelProcessorx/UpdateCounters.cpp\
    … …  
    249250                                        ./modules/ModelProcessorx/CreateParameters.cpp\
    250251                                        ./modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp\
     252                                        ./modules/ModelProcessorx/CreateFaces.cpp\
    251253                                        ./modules/ModelProcessorx/CreateEdges.cpp\
    252                                         ./modules/ModelProcessorx/CreateElementToEdgeConnectivity.cpp\
    253254                                        ./modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp\
    254255                                        ./modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp\

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r15595 r15611  
    5050
    5151                /*intialize inputs and results: */
    52                 this->inputs=new Inputs();
    53                 this->results=new Results();
     52                this->inputs  = new Inputs();
     53                this->results = new Results();
    5454
    5555                /*initialize pointers:*/

  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r15461 r15611  
    3333        this->numberofvertices=-1;
    3434        this->numberofelements=-1;
     35        this->numberoffaces=-1;
    3536        this->numberofedges=-1;
    3637        this->elements=NULL;
     38        this->faces=NULL;
    3739        this->edges=NULL;
    3840        this->elementtoedgeconnectivity      =NULL;
    … …  
    7577        FetchData(&this->numberofelements,MeshNumberofelementsEnum);
    7678        FetchData(&this->elements,NULL,NULL,MeshElementsEnum);
     79        this->faces                           = NULL;
    7780        this->edges                           = NULL;
    7881        this->elementtoedgeconnectivity       = NULL;
    … …  
    109112
    110113        xDelete<int>(this->elements);
     114        xDelete<int>(this->faces);
    111115        xDelete<int>(this->edges);
    112116        xDelete<int>(this->elementtoedgeconnectivity);

  • issm/trunk-jpl/src/c/classes/IoModel.h

    r15461 r15611  
    3434                int   numberofvertices;
    3535                int   numberofelements;
     36                int   numberoffaces;
    3637                int   numberofedges;
    3738                int  *elements;
     39                int  *faces;
    3840                int  *edges;
    3941                int  *elementtoedgeconnectivity;

  • issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp

    r15430 r15611  
    5252
    5353        /*Get edge*/
    54         int i1 = iomodel->edges[4*index+0];
    55         int i2 = iomodel->edges[4*index+1];
    56         int e1 = iomodel->edges[4*index+2];
    57         int e2 = iomodel->edges[4*index+3];
     54        int i1 = iomodel->faces[4*index+0];
     55        int i2 = iomodel->faces[4*index+1];
     56        int e1 = iomodel->faces[4*index+2];
     57        int e2 = iomodel->faces[4*index+3];
    5858
    5959        /*First, see wether this is an internal or boundary edge (if e2=-1)*/

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

    r15489 r15611  
    5757                        for(i=0;i<iomodel->numberofedges;i++){
    5858                                if(my_edges[i]){
    59                                         v1 = iomodel->edges[4*i+0]-1;
    60                                         v2 = iomodel->edges[4*i+1]-1;
     59                                        v1 = iomodel->edges[2*i+0]-1;
     60                                        v2 = iomodel->edges[2*i+1]-1;
    6161                                        if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    6262                                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
    … …  
    9999                        for(i=0;i<iomodel->numberofedges;i++){
    100100                                if(my_edges[i]){
    101                                         v1 = iomodel->edges[4*i+0]-1;
    102                                         v2 = iomodel->edges[4*i+1]-1;
     101                                        v1 = iomodel->edges[2*i+0]-1;
     102                                        v2 = iomodel->edges[2*i+1]-1;
    103103                                        values=xNew<IssmDouble>(N);
    104104                                        spcpresent=false;

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp

    r15465 r15611  
    2222        if (stabilization==3){
    2323
    24                 /*Get edges and elements*/
    25                 CreateEdges(iomodel);
     24                /*Get faces and elements*/
     25                CreateFaces(iomodel);
    2626                iomodel->FetchData(1,ThicknessEnum);
    2727
    2828                /*First load data:*/
    29                 for(int i=0;i<iomodel->numberofedges;i++){
     29                for(int i=0;i<iomodel->numberoffaces;i++){
    3030
    3131                        /*Get left and right elements*/
    32                         element=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
     32                        element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
    3333
    3434                        /*Now, if this element is not in the partition, pass: */

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateEdges.cpp

    r15435 r15611  
    1212
    1313        /*Check Iomodel properties*/
    14         if(iomodel->dim!=2)             _error_("only 2d model are supported");
    1514        if(iomodel->numberofvertices<3) _error_("not enough elements in mesh");
    1615        _assert_(iomodel->elements);
    … …  
    1918        bool exist;
    2019        int  i,j,v1,v2,v3;
    21         int  maxnbe,nbe;
     20        int  maxnbe,nbe,elementnbe,elementnbv;
     21        int *elementedges = NULL;
     22
     23        /*Mesh dependent variables*/
     24        if(iomodel->dim==2){
     25                elementnbv = 3;
     26                elementnbe = 3;
     27                elementedges=xNew<int>(elementnbe*2);
     28                elementedges[2*0+0] = 1; elementedges[2*0+1] = 2;
     29                elementedges[2*1+0] = 2; elementedges[2*1+1] = 0;
     30                elementedges[2*2+0] = 0; elementedges[2*2+1] = 1;
     31        }
     32        else if(iomodel->dim==3){
     33                elementnbv = 6;
     34                elementnbe = 9;
     35                elementedges=xNew<int>(elementnbe*2);
     36                elementedges[2*0+0] = 0; elementedges[2*0+1] = 3;
     37                elementedges[2*1+0] = 1; elementedges[2*1+1] = 4;
     38                elementedges[2*2+0] = 2; elementedges[2*2+1] = 5;
     39                elementedges[2*3+0] = 1; elementedges[2*3+1] = 2;
     40                elementedges[2*4+0] = 2; elementedges[2*4+1] = 0;
     41                elementedges[2*5+0] = 0; elementedges[2*5+1] = 1;
     42                elementedges[2*6+0] = 4; elementedges[2*6+1] = 5;
     43                elementedges[2*7+0] = 5; elementedges[2*7+1] = 3;
     44                elementedges[2*8+0] = 3; elementedges[2*8+1] = 4;
     45        }
     46        else{
     47                _error_("mesh dimension not supported yet");
     48        }
    2249
    2350        /*Maximum number of edges*/
    24         maxnbe = 3*iomodel->numberofelements;
     51        maxnbe = elementnbe*iomodel->numberofelements;
    2552
    2653        /*Initialize intermediaries*/
    27         int*  edgestemp = xNew<int>(maxnbe*4);         /*format: [vertex1 vertex2 element1 element2]                */
    28         bool* exchange  = xNewZeroInit<bool>(maxnbe);  /*Edges are ordered, we need to keep track of vertex swapping*/
    29         for(i=0;i<maxnbe;i++) edgestemp[i*4+3]=-1;     /*Initialize last column of edges as -1 (boundary edge)      */
     54        int *edgestemp                 = xNew<int>(maxnbe*2);                             /*format: [vertex1 vertex2]       */
     55        int *element_edge_connectivity = xNew<int>(iomodel->numberofelements*elementnbe); /*format: [edge1 edge2 ... edgen] */
    3056
    3157        /*Initialize chain*/
    … …  
    3864
    3965        for(i=0;i<iomodel->numberofelements;i++){
    40                 for(j=0;j<3;j++){
     66                for(j=0;j<elementnbe;j++){
    4167
    42                         /*Get the two indices of the edge number j of the ith triangle*/
    43                         v1 = iomodel->elements[i*3+j];
    44                         if(j==2)
    45                          v2 = iomodel->elements[i*3+0];
    46                         else
    47                          v2 = iomodel->elements[i*3+j+1];
     68                        /*Get the two indices of the edge number j of the ith element*/
     69                        v1 = iomodel->elements[i*elementnbv+elementedges[2*j+0]]-1; _assert_(v1>=0 & v1<iomodel->numberofvertices);
     70                        v2 = iomodel->elements[i*elementnbv+elementedges[2*j+1]]-1; _assert_(v2>=0 & v2<iomodel->numberofvertices);
    4871
    4972                        /*v1 and v2 must be sorted*/
    … …  
    5679
    5780                        /*Go through all processed edges connected to v1 and check whether we have seen this edge yet*/
    58                         _assert_(v1>=0 & v1<iomodel->numberofvertices);
    5981                        for(int e=head_minv[v1]; e!=-1; e=next_edge[e]){
    60                                 if(edgestemp[e*4+1]==v2){
     82                                if(edgestemp[e*2+1]==v2+1){
    6183                                        exist = true;
    62                                         edgestemp[e*4+3]=i+1;
     84                                        element_edge_connectivity[i*elementnbv+j]=e;
    6385                                        break;
    6486                                }
    … …  
    7092
    7193                                /*Update edges*/
    72                                 edgestemp[nbe*4+0] = v1;
    73                                 edgestemp[nbe*4+1] = v2;
    74                                 edgestemp[nbe*4+2] = i+1;
    75                                 if(v1!=iomodel->elements[i*3+j]) exchange[nbe]=true;
     94                                edgestemp[nbe*2+0] = v1+1;
     95                                edgestemp[nbe*2+1] = v2+1;
     96
     97                                /*Update Connectivity*/
     98                                element_edge_connectivity[i*elementnbv+j]=nbe;
    7699
    77100                                /*Update chain*/
    … …  
    90113
    91114        /*Create final edges*/
    92         int* edges = xNew<int>(nbe*4); /*vertex1 vertex2 element1 element2*/
    93         for(int i=0;i<nbe;i++){
    94                 if(exchange[i]){
    95                         edges[i*4+0]=edgestemp[i*4+1];
    96                         edges[i*4+1]=edgestemp[i*4+0];
    97                 }
    98                 else{
    99                         edges[i*4+0]=edgestemp[i*4+0];
    100                         edges[i*4+1]=edgestemp[i*4+1];
    101                 }
    102                 edges[i*4+2]=edgestemp[i*4+2];
    103                 edges[i*4+3]=edgestemp[i*4+3];
    104         }
     115        int* edges = xNew<int>(nbe*2); /*vertex1 vertex2*/
     116        for(int i=0;i<2*nbe;i++) edges[i] = edgestemp[i];
     117
     118        /*Clean up*/
    105119        xDelete<int>(edgestemp);
    106         xDelete<bool>(exchange);
     120        xDelete<int>(elementedges);
    107121
    108122        /*Assign output pointers*/
    109         iomodel->edges         = edges;
    110         iomodel->numberofedges = nbe;
     123        iomodel->edges                     = edges;
     124        iomodel->elementtoedgeconnectivity = element_edge_connectivity;
     125        iomodel->numberofedges             = nbe;
    111126}

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp

    r15583 r15611  
    4545                case P2Enum:
    4646                        EdgesPartitioning(&my_edges,iomodel);
    47                         CreateElementToEdgeConnectivity(iomodel);
    4847                        for(i=0;i<iomodel->numberofvertices;i++){
    4948                                if(iomodel->my_vertices[i]){

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

    r15567 r15611  
    325325                        if(my_edges[i]){
    326326
    327                                 v1 = iomodel->edges[4*i+0]-1;
    328                                 v2 = iomodel->edges[4*i+1]-1;
     327                                v1 = iomodel->edges[2*i+0]-1;
     328                                v2 = iomodel->edges[2*i+1]-1;
    329329
    330330                                if(!xIsNan<IssmDouble>(spcvx[v1]) && !xIsNan<IssmDouble>(spcvx[v2])){

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/EdgesPartitioning.cpp

    r15461 r15611  
    1111
    1212        /*Intermediaries*/
    13         int  el1,el2;
    14         bool my_edge;
     13        int elementnbe;
    1514
    1615        /*Get edges and elements*/
    1716        CreateEdges(iomodel);
     17        _assert_(iomodel->elementtoedgeconnectivity);
    1818
     19        /*Mesh dependent variables*/
     20        if(iomodel->dim==2){
     21                elementnbe = 3;
     22        }
     23        else if(iomodel->dim==3){
     24                elementnbe = 9;
     25        }
     26        else{
     27                _error_("mesh dimension not supported yet");
     28        }
    1929        /*output: */
    20         bool* my_edges=xNew<bool>(iomodel->numberofedges);
     30        bool* my_edges=xNewZeroInit<bool>(iomodel->numberofedges);
    2131
    22         for(int i=0;i<iomodel->numberofedges;i++){
    23 
    24                 /*Get left and right elements*/
    25                 el1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
    26                 el2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
    27 
    28                 /*Check whether we should include this edge (el2 is -2 for boundary edges)*/
    29                 my_edge = iomodel->my_elements[el1];
    30                 if(!my_edge && el2>=0){
    31                         my_edge = iomodel->my_elements[el2];
     32        for(int i=0;i<iomodel->numberofelements;i++){
     33                if(iomodel->my_elements[i]){
     34                        for(int j=0;j<elementnbe;j++){
     35                                my_edges[iomodel->elementtoedgeconnectivity[i*elementnbe+j]] = true;
     36                        }
    3237                }
    33 
    34                 my_edges[i] = my_edge;
    3538        }
    3639

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r15583 r15611  
    117117void ElementsAndVerticesPartitioning(bool** pmy_elements, int** pmy_vertices, IoModel* iomodel);
    118118void NodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices,  IoModel* iomodel, bool continuous);
    119 void EdgesPartitioning(bool** pmy_nodes,IoModel* iomodel);
     119void FacesPartitioning(bool** pmy_faces,IoModel* iomodel);
     120void EdgesPartitioning(bool** pmy_edges,IoModel* iomodel);
     121
     122/*Mesh properties*/
     123void CreateFaces(IoModel* iomodel);
     124void CreateEdges(IoModel* iomodel);
    120125
    121126/*Connectivity*/
    122 void CreateEdges(IoModel* iomodel);
    123 void CreateElementToEdgeConnectivity(IoModel* iomodel);
    124127void CreateSingleNodeToElementConnectivity(IoModel* iomodel);
    125128void CreateNumberNodeToElementConnectivity(IoModel* iomodel);

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp

    r15461 r15611  
    3939
    4040        /* Each element has it own nodes (as many as vertices) + additional nodes
    41          * from neighboring elements for each edge. This yields to a very different
     41         * from neighboring elements for each face. This yields to a very different
    4242         * partition for the nodes and the vertices. The vertices are similar to
    43          * continuous galerkin, but the nodes partitioning involves edges, which
    44          * mess up sorting of ids. */
     43         * continuous galerkin, but the nodes partitioning involves faces, which
     44         * messes up sorting of ids. */
    4545
    4646
    … …  
    5151        int  pos;
    5252
    53         /*Get edges and elements*/
     53        /*Get faces and elements*/
    5454        CreateEdges(iomodel);
    5555
    … …  
    5757         *  - there are three nodes per element (discontinous)
    5858         *  - for each element present of each partition, its three nodes will be in this partition
    59          *  - the edges require the dofs of the 2 nodes of each elements sharing the edge.
    60          *    if the 2 elements sharing the edge are on 2 different cpus, we must duplicate
    61          *    the two nodes that are not on the cpus so that the edge can access the dofs of
     59         *  - the faces require the dofs of the 2 nodes of each elements sharing the face.
     60         *    if the 2 elements sharing the face are on 2 different cpus, we must duplicate
     61         *    the two nodes that are not on the cpus so that the face can access the dofs of
    6262         *    all its 4 nodes
    6363         */
    … …  
    6767
    6868        /*First: add all the nodes of all the elements belonging to this cpu*/
    69         if (iomodel->dim==2){
     69        if(iomodel->dim==2){
    7070                for (i=0;i<iomodel->numberofelements;i++){
    7171                        if (my_elements[i]){
    … …  
    8282        /*Second: add all missing nodes*/
    8383
    84         /*Get edges and elements*/
    85         CreateEdges(iomodel);
     84        /*Get faces and elements*/
     85        CreateFaces(iomodel);
    8686
    8787        /*!All elements have been partitioned above, only create elements for this CPU: */
    88         for(int i=0;i<iomodel->numberofedges;i++){
     88        for(int i=0;i<iomodel->numberoffaces;i++){
    8989
    9090                /*Get left and right elements*/
    91                 e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
    92                 e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
     91                e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     92                e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
    9393
    9494                /* 1) If the element e1 is in the current partition
    95                  * 2) and if the edge of the element is shared by another element (internal edge)
     95                 * 2) and if the face of the element is shared by another element (internal face)
    9696                 * 3) and if this element is not in the same partition:
    9797                 * we must clone the nodes on this partition so that the loads (Numericalflux)
    … …  
    100100
    101101                        /*1: Get vertices ids*/
    102                         i1=iomodel->edges[4*i+0];
    103                         i2=iomodel->edges[4*i+1];
     102                        i1=iomodel->faces[4*i+0];
     103                        i2=iomodel->faces[4*i+1];
    104104
    105105                        /*2: Get the column where these ids are located in the index*/
    … …  
    124124                        }
    125125                        else{
    126                                 _error_("Problem in edges creation");
     126                                _error_("Problem in faces creation");
    127127                        }
    128128                }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp

    r15465 r15611  
    2525        if (stabilization==3){
    2626
    27                 /*Get edges and elements*/
    28                 CreateEdges(iomodel);
     27                /*Get faces and elements*/
     28                CreateFaces(iomodel);
    2929                iomodel->FetchData(1,ThicknessEnum);
    3030
    3131                /*First load data:*/
    32                 for(int i=0;i<iomodel->numberofedges;i++){
     32                for(int i=0;i<iomodel->numberoffaces;i++){
    3333
    3434                        /*Get left and right elements*/
    35                         element=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
     35                        element=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
    3636
    3737                        /*Now, if this element is not in the partition, pass: */
Note: See TracChangeset for help on using the changeset viewer.