Changeset 23576


Ignore:
Timestamp:
12/28/18 09:17:00 (6 years ago)
Author:
Mathieu Morlighem
Message:

CHG: changing way clones are flagged and nodes are numbered

Location:
issm/trunk-jpl/src/c
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r23561 r23576  
    545545        else{
    546546                /*Coupling: we are going to create P1 Elements only*/
    547 
    548                 Node*  node  = NULL;
    549                 int    lid=0;
    550                 if(!nodes) nodes = new Nodes();
    551 
    552547                iomodel->FetchData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS",
    553548                                        "md.flowequation.vertex_equation","md.stressbalance.referential");
    554549                if(isFS){
    555                         /*P1+ velocity*/
     550                        int* approximations = xNew<int>(2*iomodel->numberofvertices+iomodel->numberofelements);
    556551                        for(int i=0;i<iomodel->numberofvertices;i++){
    557                                 if(iomodel->my_vertices[i]){
    558                                         approximation=IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]));
    559                                         if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
    560                                         nodes->AddObject(new Node(i+1,i,lid++,0,i,false,iomodel,StressbalanceAnalysisEnum,approximation));
    561                                 }
     552                                approximation = IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]));
     553                                if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
     554                                approximations[i] = approximation;
    562555                        }
    563                         for(int i=0;i<iomodel->numberofelements;i++){
    564                                 if(iomodel->my_elements[i]){
    565                                         node = new Node(iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
    566                                         node->Deactivate();
    567                                         nodes->AddObject(node);
    568                                 }
    569                         }
    570                         /*P1 pressure*/
    571                         for(int i=0;i<iomodel->numberofvertices;i++){
    572                                 if(iomodel->my_vertices[i]){
    573                                         approximation=IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]));
    574                                         node = new Node(iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,0,i,false,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
     556                        for(int i=0;i<iomodel->numberofelements;i++) approximations[iomodel->numberofvertices+i] = FSvelocityEnum;
     557                        for(int i=0;i<iomodel->numberofvertices;i++) approximations[iomodel->numberofvertices+iomodel->numberofelements+i] = FSpressureEnum;
     558                        ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,MINIcondensedEnum,0,approximations);
     559                        xDelete<int>(approximations);
     560
     561                        for(int i=0;i<nodes->Size();i++){
     562                                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
     563                                int   sid = node->Sid();
     564                                if(sid>=iomodel->numberofvertices+iomodel->numberofelements){
     565                                        /*Constrain pressure if not FS*/
     566                                        int id = sid - (iomodel->numberofvertices+iomodel->numberofelements);
     567                                        approximation=IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[id]));
    575568                                        if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
    576569                                                node->Deactivate();
    577570                                        }
    578                                         nodes->AddObject(node);
    579571                                }
    580572                        }
    581573                }
    582574                else{
    583                         for(int i=0;i<iomodel->numberofvertices;i++){
    584                                 if(iomodel->my_vertices[i]){
    585                                         nodes->AddObject(new Node(i+1,i,lid++,0,i,false,iomodel,StressbalanceAnalysisEnum,IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]))));
    586                                 }
    587                         }
     575                        int* approximations = xNew<int>(iomodel->numberofvertices);
     576                        for(int i=0;i<iomodel->numberofvertices;i++) approximations[i] = IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]));
     577                        ::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,P1Enum,0,approximations);
     578                        xDelete<int>(approximations);
    588579                }
    589580                iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS",

  • issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp

    r23561 r23576  
    8080        /*Intermediaries*/
    8181        bool  isSIA;
    82         Node* node = NULL;
    8382
    8483        /*Fetch parameters: */
    … …  
    9594        }
    9695
    97         for(int i=0;i<iomodel->numberofvertices;i++){
    98                 if(iomodel->my_vertices[i]){
    99 
    100                         /*Create new node if is in this processor's partition*/
    101                         node = new Node(i+1,i,lid++,0,i,false,iomodel,StressbalanceSIAAnalysisEnum,IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i])));
    102 
    103                         /*Deactivate node if not SIA*/
    104                         if(IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[i]))!=SIAApproximationEnum){
    105                                 node->Deactivate();
    106                         }
    107 
    108                         /*Add to Nodes dataset*/
    109                         nodes->AddObject(node);
    110                 }
    111         }
    112 
     96        ::CreateNodes(nodes,iomodel,StressbalanceSIAAnalysisEnum,P1Enum,SIAApproximationEnum);
     97        for(int i=0;i<nodes->Size();i++){
     98                Node* node=xDynamicCast<Node*>(nodes->GetObjectByOffset(i));
     99                int   sid = node->Sid();
     100                int approximation=IoCodeToEnumVertexEquation(reCast<int>(iomodel->Data("md.flowequation.vertex_equation")[sid]));
     101                if(approximation!=SIAApproximationEnum) node->Deactivate();
     102        }
    113103        iomodel->DeleteData(6,"md.mesh.vertexonbase","md.mesh.vertexonsurface","md.flowequation.borderSSA","md.flowequation.borderFS","md.flowequation.vertex_equation","md.stressbalance.referential");
    114104

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

    r23532 r23576  
    29182918                        penta_node_ids[ 4]=iomodel->elements[6*index+4];
    29192919                        penta_node_ids[ 5]=iomodel->elements[6*index+5];
    2920                         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
    2921                         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
    2922                         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
     2920                        penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+0]+1;
     2921                        penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+1]+1;
     2922                        penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*index+2]+1;
    29232923                        break;
    29242924                case P1xP3Enum:
    … …  
    29312931                        penta_node_ids[ 4]=iomodel->elements[6*index+4];
    29322932                        penta_node_ids[ 5]=iomodel->elements[6*index+5];
    2933                         penta_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+0]+1;
    2934                         penta_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+1]+1;
    2935                         penta_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+2]+1;
    2936                         penta_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+0]+2;
    2937                         penta_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+1]+2;
    2938                         penta_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoedgeconnectivity[9*index+2]+2;
     2933                        penta_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1;
     2934                        penta_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1;
     2935                        penta_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1;
     2936                        penta_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2;
     2937                        penta_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2;
     2938                        penta_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2;
    29392939                        break;
    29402940                case P2xP1Enum:
    … …  
    29472947                        penta_node_ids[ 4]=iomodel->elements[6*index+4];
    29482948                        penta_node_ids[ 5]=iomodel->elements[6*index+5];
    2949                         penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
    2950                         penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
    2951                         penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
    2952                         penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
    2953                         penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
    2954                         penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
     2949                        penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+0]+1;
     2950                        penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+1]+1;
     2951                        penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+2]+1;
     2952                        penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+3]+1;
     2953                        penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+4]+1;
     2954                        penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*index+5]+1;
    29552955                        break;
    29562956                case P1xP4Enum:
    … …  
    29632963                        penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/
    29642964                        penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/
    2965                         penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
    2966                         penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
    2967                         penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
    2968                         penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/
    2969                         penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/
    2970                         penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/
    2971                         penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 3/4 vertical edge 1*/
    2972                         penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 3/4 vertical edge 2*/
    2973                         penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 3/4 vertical edge 3*/
     2965                        penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /*mid vertical edge 1*/
     2966                        penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /*mid vertical edge 2*/
     2967                        penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /*mid vertical edge 3*/
     2968                        penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 1/4 vertical edge 1*/
     2969                        penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 1/4 vertical edge 2*/
     2970                        penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 1/4 vertical edge 3*/
     2971                        penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+0]+3; /* 3/4 vertical edge 1*/
     2972                        penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+1]+3; /* 3/4 vertical edge 2*/
     2973                        penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*index+2]+3; /* 3/4 vertical edge 3*/
    29742974                        break;
    29752975                case P2xP4Enum:
    … …  
    29822982                        penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/
    29832983                        penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/
    2984                         penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
    2985                         penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
    2986                         penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
    2987                         penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
    2988                         penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
    2989                         penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
    2990                         penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
    2991                         penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
    2992                         penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
    2993                         penta_node_ids[15]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/
    2994                         penta_node_ids[16]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/
    2995                         penta_node_ids[17]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/
    2996                         penta_node_ids[18]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 3/4 vertical edge 1*/
    2997                         penta_node_ids[19]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 3/4 vertical edge 2*/
    2998                         penta_node_ids[20]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 3/4 vertical edge 3*/
    2999                         penta_node_ids[21]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+1; /* 1/4 vertical face 1*/
    3000                         penta_node_ids[22]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+1; /* 1/4 vertical face 2*/
    3001                         penta_node_ids[23]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+1; /* 1/4 vertical face 3*/
    3002                         penta_node_ids[24]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+2; /* 2/4 vertical face 1*/
    3003                         penta_node_ids[25]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+2; /* 2/4 vertical face 2*/
    3004                         penta_node_ids[26]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+2; /* 2/4 vertical face 3*/
    3005                         penta_node_ids[27]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+3; /* 3/4 vertical face 1*/
    3006                         penta_node_ids[28]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+3; /* 3/4 vertical face 2*/
    3007                         penta_node_ids[29]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+3; /* 3/4 vertical face 3*/
     2984                        penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
     2985                        penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
     2986                        penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
     2987                        penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
     2988                        penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
     2989                        penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
     2990                        penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
     2991                        penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
     2992                        penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
     2993                        penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /* 1/4 vertical edge 1*/
     2994                        penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /* 1/4 vertical edge 2*/
     2995                        penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /* 1/4 vertical edge 3*/
     2996                        penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 3/4 vertical edge 1*/
     2997                        penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 3/4 vertical edge 2*/
     2998                        penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 3/4 vertical edge 3*/
     2999                        penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; /* 1/4 vertical face 1*/
     3000                        penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; /* 1/4 vertical face 2*/
     3001                        penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; /* 1/4 vertical face 3*/
     3002                        penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+2; /* 2/4 vertical face 1*/
     3003                        penta_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+2; /* 2/4 vertical face 2*/
     3004                        penta_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+2; /* 2/4 vertical face 3*/
     3005                        penta_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+3; /* 3/4 vertical face 1*/
     3006                        penta_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+3; /* 3/4 vertical face 2*/
     3007                        penta_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+3; /* 3/4 vertical face 3*/
    30083008                        break;
    30093009                case P2Enum:
    … …  
    30253025                        penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
    30263026                        penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
    3027                         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
    3028                         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
    3029                         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
     3027                        penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+0]+1;
     3028                        penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+1]+1;
     3029                        penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+2]+1;
    30303030                        break;
    30313031                case P2bubbleEnum: case P2bubblecondensedEnum:
    … …  
    31053105                        penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
    31063106                        penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
    3107                         penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
    3108                         penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
    3109                         penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
    3110 
    3111                         penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+0];
    3112                         penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+1];
    3113                         penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+2];
    3114                         penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+3];
    3115                         penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+4];
    3116                         penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->elements[6*index+5];
     3107                        penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+0]+1;
     3108                        penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+1]+1;
     3109                        penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*index+2]+1;
     3110
     3111                        penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+0];
     3112                        penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+1];
     3113                        penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+2];
     3114                        penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+3];
     3115                        penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+4];
     3116                        penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*index+5];
    31173117                        break;
    31183118                case LATaylorHoodEnum:
    … …  
    31473147                        penta_node_ids[ 4]=iomodel->elements[6*index+4]; /*Vertex 5*/
    31483148                        penta_node_ids[ 5]=iomodel->elements[6*index+5]; /*Vertex 6*/
    3149                         penta_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
    3150                         penta_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
    3151                         penta_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
    3152                         penta_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
    3153                         penta_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
    3154                         penta_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
    3155                         penta_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
    3156                         penta_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
    3157                         penta_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
    3158                         penta_node_ids[15]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+2; /* 1/4 vertical edge 1*/
    3159                         penta_node_ids[16]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+2; /* 1/4 vertical edge 2*/
    3160                         penta_node_ids[17]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+2; /* 1/4 vertical edge 3*/
    3161                         penta_node_ids[18]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+0]+3; /* 2/4 vertical edge 1*/
    3162                         penta_node_ids[19]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+1]+3; /* 2/4 vertical edge 2*/
    3163                         penta_node_ids[20]=iomodel->numberofvertices+3*iomodel->elementtoedgeconnectivity[9*index+2]+3; /* 2/4 vertical edge 3*/
    3164                         penta_node_ids[21]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+1; /* 1/4 vertical face 1*/
    3165                         penta_node_ids[22]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+1; /* 1/4 vertical face 2*/
    3166                         penta_node_ids[23]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+1; /* 1/4 vertical face 3*/
    3167                         penta_node_ids[24]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+2; /* 2/4 vertical face 1*/
    3168                         penta_node_ids[25]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+2; /* 2/4 vertical face 2*/
    3169                         penta_node_ids[26]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+2; /* 2/4 vertical face 3*/
    3170                         penta_node_ids[27]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+2]+3; /* 3/4 vertical face 1*/
    3171                         penta_node_ids[28]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+3]+3; /* 3/4 vertical face 2*/
    3172                         penta_node_ids[29]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->elementtofaceconnectivity[5*index+4]+3; /* 3/4 vertical face 3*/
    3173 
    3174                         penta_node_ids[30]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+0];
    3175                         penta_node_ids[31]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+1];
    3176                         penta_node_ids[32]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+2];
    3177                         penta_node_ids[33]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+3];
    3178                         penta_node_ids[34]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+4];
    3179                         penta_node_ids[35]=iomodel->numberofvertices+3*iomodel->numberofedges+3*iomodel->numberoffaces+iomodel->elements[6*index+5];
     3149                        penta_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1; /*mid vertical edge 1*/
     3150                        penta_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1; /*mid vertical edge 2*/
     3151                        penta_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1; /*mid vertical edge 3*/
     3152                        penta_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1; /*mid basal edge 1*/
     3153                        penta_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1; /*mid basal edge 2*/
     3154                        penta_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1; /*mid basal edge 3*/
     3155                        penta_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1; /*mid top edge 1*/
     3156                        penta_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1; /*mid top edge 2*/
     3157                        penta_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; /*mid top edge 3*/
     3158                        penta_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+1; /* 1/4 vertical edge 1*/
     3159                        penta_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+1; /* 1/4 vertical edge 2*/
     3160                        penta_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+1; /* 1/4 vertical edge 3*/
     3161                        penta_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+0]+2; /* 3/4 vertical edge 1*/
     3162                        penta_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+1]+2; /* 3/4 vertical edge 2*/
     3163                        penta_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*index+2]+2; /* 3/4 vertical edge 3*/
     3164                        penta_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+1; /* 1/4 vertical face 1*/
     3165                        penta_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+1; /* 1/4 vertical face 2*/
     3166                        penta_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+1; /* 1/4 vertical face 3*/
     3167                        penta_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+2; /* 2/4 vertical face 1*/
     3168                        penta_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+2; /* 2/4 vertical face 2*/
     3169                        penta_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+2; /* 2/4 vertical face 3*/
     3170                        penta_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+0]+3; /* 3/4 vertical face 1*/
     3171                        penta_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+1]+3; /* 3/4 vertical face 2*/
     3172                        penta_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*index+2]+3; /* 3/4 vertical face 3*/
     3173
     3174                        penta_node_ids[30]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+0];
     3175                        penta_node_ids[31]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+1];
     3176                        penta_node_ids[32]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+2];
     3177                        penta_node_ids[33]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+3];
     3178                        penta_node_ids[34]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+4];
     3179                        penta_node_ids[35]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*index+5];
    31803180                        break;
    31813181                case CrouzeixRaviartEnum:

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

    r23561 r23576  
    27072707                }
    27082708                SpcNodesx(new_nodes_list[i],new_constraints_list[i],this->parameters,analysis_type);
     2709                new_nodes_list[i]->FlagClones(analysis_type);/*FIXME: should be removed !*/
    27092710                NodesDofx(new_nodes_list[i],this->parameters,analysis_type);
    27102711        }

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

    r23532 r23576  
    113113        }
    114114
    115         int count=0;
     115        int count =0;
    116116        if(M==iomodel->numberofvertices){
    117117                switch(finite_element){
    … …  
    127127                                break;
    128128                        case P2Enum:
     129                                for(i=0;i<iomodel->numberofvertices;i++){
     130                                        if((iomodel->my_vertices[i])){
     131                                                if (!xIsNan<IssmDouble>(spcdata[i])){
     132                                                        constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
     133                                                        count++;
     134                                                }
     135                                        }
     136                                }
     137                                for(i=0;i<iomodel->numberofedges;i++){
     138                                        if(iomodel->my_edges[i] && boundaryedge[i]){
     139                                                v1 = iomodel->edges[3*i+0]-1;
     140                                                v2 = iomodel->edges[3*i+1]-1;
     141                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     142                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
     143                                                                                        dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
     144                                                        count++;
     145                                                }
     146                                        }
     147                                }
     148                                if(iomodel->meshelementtype==PentaEnum){
     149                                        for(i=0;i<iomodel->numberofverticalfaces;i++){
     150                                                if(iomodel->my_vfaces[i]){
     151                                                        value=0.;
     152                                                        for(j=0;j<4;j++) value += spcdata[iomodel->verticalfaces[i*4+j] -1]/4.;
     153                                                        if(!xIsNan<IssmDouble>(value)){
     154                                                                constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+iomodel->numberofedges+i+1,dof,value,analysis_type));
     155                                                                count++;
     156                                                        }
     157                                                }
     158                                        }
     159                                }
     160                                break;
     161                        case P2bubbleEnum:
    129162                                for(i=0;i<iomodel->numberofvertices;i++){
    130163                                        if((iomodel->my_vertices[i])){
    … …  
    166199                                        }
    167200                                }
    168                                 break;
    169                         case P2bubbleEnum:
    170                                 for(i=0;i<iomodel->numberofvertices;i++){
    171                                         if((iomodel->my_vertices[i])){
    172                                                 if (!xIsNan<IssmDouble>(spcdata[i])){
    173                                                         constraints->AddObject(new SpcStatic(count+1,i+1,dof,spcdata[i],analysis_type));
    174                                                         count++;
    175                                                 }
    176                                         }
    177                                 }
    178                                 for(i=0;i<iomodel->numberofedges;i++){
    179                                         if(iomodel->my_edges[i] && boundaryedge[i]){
    180                                                 v1 = iomodel->edges[3*i+0]-1;
    181                                                 v2 = iomodel->edges[3*i+1]-1;
    182                                                 if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    183                                                         constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
    184                                                                                         dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
    185                                                         count++;
    186                                                 }
    187                                         }
    188                                 }
    189                                 if(iomodel->meshelementtype==PentaEnum){
    190                                         for(i=0;i<iomodel->numberoffaces;i++){
    191                                                 if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    192                                                         if(iomodel->my_faces[i]){
    193                                                                 numfacevertices = iomodel->faces[i*iomodel->facescols+3];
    194                                                                 value=0.;
    195                                                                 for(j=0;j<numfacevertices;j++){
    196                                                                         value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1];
    197                                                                 }
    198                                                                 value = value/reCast<IssmDouble>(numfacevertices);
    199                                                                 if(!xIsNan<IssmDouble>(value)){
    200                                                                         constraints->AddObject(new SpcStatic(count+1,
    201                                                                                                         iomodel->numberofvertices+iomodel->numberofedges+i+1,
    202                                                                                                         dof,value,analysis_type));
    203                                                                         count++;
    204                                                                 }
    205                                                         }
    206                                                 }
    207                                         }
    208                                 }
    209201                                for(i=0;i<iomodel->numberofelements;i++){
    210202                                        if(iomodel->my_elements[i]){
    … …  
    233225                                }
    234226                                for(i=0;i<iomodel->numberofedges;i++){
    235                                         if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
    236                                                 if(iomodel->my_edges[i]){
    237                                                         v1 = iomodel->edges[3*i+0]-1;
    238                                                         v2 = iomodel->edges[3*i+1]-1;
    239                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    240                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+3*i+1,
    241                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    242                                                                 constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2,
    243                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    244                                                                 constraints->AddObject(new SpcStatic(count+3,iomodel->numberofvertices+3*i+3,
    245                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    246                                                                 count=count+3;
    247                                                         }
    248                                                 }
    249                                         }
    250                                         if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
    251                                                 if(iomodel->my_edges[i]){
    252                                                         v1 = iomodel->edges[3*i+0]-1;
    253                                                         v2 = iomodel->edges[3*i+1]-1;
    254                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    255                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+3*i+1,
    256                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    257                                                                 count=count+1;
    258                                                         }
    259                                                 }
    260                                         }
    261                                 }
    262                                 for(i=0;i<iomodel->numberoffaces;i++){
    263                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    264                                                 if(iomodel->my_faces[i]){
    265                                                         numfacevertices = iomodel->faces[i*iomodel->facescols+3];
    266                                                         value=0.;
    267                                                         for(j=0;j<numfacevertices;j++){
    268                                                                 value += spcdata[iomodel->faces[i*iomodel->facescols+4+j] -1];
    269                                                         }
    270                                                         value = value/reCast<IssmDouble>(numfacevertices);
    271                                                         if(!xIsNan<IssmDouble>(value)){
    272                                                                 constraints->AddObject(new SpcStatic(count+1,
    273                                                                                                 iomodel->numberofvertices+3*iomodel->numberofedges+3*i+1,
    274                                                                                                 dof,value,analysis_type));
    275                                                                 constraints->AddObject(new SpcStatic(count+2,
    276                                                                                                 iomodel->numberofvertices+3*iomodel->numberofedges+3*i+2,
    277                                                                                                 dof,value,analysis_type));
    278                                                                 constraints->AddObject(new SpcStatic(count+3,
    279                                                                                                 iomodel->numberofvertices+3*iomodel->numberofedges+3*i+3,
    280                                                                                                 dof,value,analysis_type));
    281                                                                 count=count+3;
    282                                                         }
     227                                        if(iomodel->my_edges[i]){
     228                                                v1 = iomodel->edges[3*i+0]-1;
     229                                                v2 = iomodel->edges[3*i+1]-1;
     230                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     231                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
     232                                                                                        dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
     233                                                        count++;
     234                                                }
     235                                        }
     236                                }
     237                                for(i=0;i<iomodel->numberofverticaledges;i++){
     238                                        if(iomodel->my_vedges[i]){
     239                                                v1 = iomodel->verticaledges[2*i+0]-1;
     240                                                v2 = iomodel->verticaledges[2*i+1]-1;
     241                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     242                                                        /*FIXME: should not be 1/2 but 1/4 and 3/4!*/
     243                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+iomodel->numberofedges+2*i+1,
     244                                                                                        dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
     245                                                        constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+iomodel->numberofedges+2*i+2,
     246                                                                                        dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
     247                                                        count=count+2;
     248                                                }
     249                                        }
     250                                }
     251                                for(i=0;i<iomodel->numberofverticalfaces;i++){
     252                                        if(iomodel->my_vfaces[i]){
     253                                                value=0.;
     254                                                for(j=0;j<4;j++) value += spcdata[iomodel->verticalfaces[i*4+j]-1]/4.;
     255                                                if(!xIsNan<IssmDouble>(value)){
     256                                                        constraints->AddObject(new SpcStatic(count+1,
     257                                                                                        iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+1,
     258                                                                                        dof,value,analysis_type));
     259                                                        constraints->AddObject(new SpcStatic(count+2,
     260                                                                                        iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+2,
     261                                                                                        dof,value,analysis_type));
     262                                                        constraints->AddObject(new SpcStatic(count+3,
     263                                                                                        iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*i+3,
     264                                                                                        dof,value,analysis_type));
     265                                                        count=count+3;
    283266                                                }
    284267                                        }
    … …  
    326309                                        }
    327310                                }
    328                                 for(i=0;i<iomodel->numberofedges;i++){
    329                                         if(iomodel->edges[i*3+2]==2){
    330                                                 if(iomodel->my_edges[i]){
    331                                                         v1 = iomodel->edges[3*i+0]-1;
    332                                                         v2 = iomodel->edges[3*i+1]-1;
    333                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    334                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
    335                                                                                                 dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
    336                                                                 count++;
    337                                                         }
     311                                for(i=0;i<iomodel->numberofverticaledges;i++){
     312                                        if(iomodel->my_vedges[i]){
     313                                                v1 = iomodel->verticaledges[2*i+0]-1;
     314                                                v2 = iomodel->verticaledges[2*i+1]-1;
     315                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     316                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
     317                                                        count++;
    338318                                                }
    339319                                        }
    … …  
    349329                                        }
    350330                                }
    351                                 for(i=0;i<iomodel->numberofedges;i++){
    352                                         if(iomodel->edges[i*3+2]==2){
    353                                                 if(iomodel->my_edges[i]){
    354                                                         v1 = iomodel->edges[3*i+0]-1;
    355                                                         v2 = iomodel->edges[3*i+1]-1;
    356                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    357                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+2*i+1,
    358                                                                                                 dof,2./3.*spcdata[v1]+1./3.*spcdata[v2],analysis_type));
    359                                                                 constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+2*i+2,
    360                                                                                                 dof,1./3.*spcdata[v1]+2./3.*spcdata[v2],analysis_type));
    361                                                                 count=count+2;
    362                                                         }
     331                                for(i=0;i<iomodel->numberofverticaledges;i++){
     332                                        if(iomodel->my_vedges[i]){
     333                                                v1 = iomodel->verticaledges[2*i+0]-1;
     334                                                v2 = iomodel->verticaledges[2*i+1]-1;
     335                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     336                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+2*i+1,
     337                                                                                        dof,2./3.*spcdata[v1]+1./3.*spcdata[v2],analysis_type));
     338                                                        constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+2*i+2,
     339                                                                                        dof,1./3.*spcdata[v1]+2./3.*spcdata[v2],analysis_type));
     340                                                        count=count+2;
    363341                                                }
    364342                                        }
    … …  
    374352                                        }
    375353                                }
    376                                 for(i=0;i<iomodel->numberofedges;i++){
    377                                         if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
    378                                                 if(iomodel->my_edges[i]){
    379                                                         v1 = iomodel->edges[3*i+0]-1;
    380                                                         v2 = iomodel->edges[3*i+1]-1;
    381                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    382                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+3*i+1,
    383                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    384                                                                 constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2,
    385                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    386                                                                 constraints->AddObject(new SpcStatic(count+3,iomodel->numberofvertices+3*i+3,
    387                                                                                                 dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
    388                                                                 count=count+3;
    389                                                         }
     354                                for(i=0;i<iomodel->numberofverticaledges;i++){
     355                                        if(iomodel->my_vedges[i]){
     356                                                v1 = iomodel->verticaledges[2*i+0]-1;
     357                                                v2 = iomodel->verticaledges[2*i+1]-1;
     358                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     359                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+3*i+1,
     360                                                                                        dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
     361                                                        constraints->AddObject(new SpcStatic(count+2,iomodel->numberofvertices+3*i+2,
     362                                                                                        dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
     363                                                        constraints->AddObject(new SpcStatic(count+3,iomodel->numberofvertices+3*i+3,
     364                                                                                        dof,1./2.*spcdata[v1]+1./2.*spcdata[v2],analysis_type));
     365                                                        count=count+3;
    390366                                                }
    391367                                        }
    … …  
    401377                                        }
    402378                                }
    403                                 for(i=0;i<iomodel->numberofedges;i++){
    404                                         if(iomodel->edges[i*3+2]!=2){
    405                                                 if(iomodel->my_edges[i]){
    406                                                         v1 = iomodel->edges[3*i+0]-1;
    407                                                         v2 = iomodel->edges[3*i+1]-1;
    408                                                         if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
    409                                                                 constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,
    410                                                                                                 dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
    411                                                                 count++;
    412                                                         }
     379                                for(i=0;i<iomodel->numberofhorizontaledges;i++){
     380                                        if(iomodel->my_hedges[i]){
     381                                                v1 = iomodel->horizontaledges[2*i+0]-1;
     382                                                v2 = iomodel->horizontaledges[2*i+1]-1;
     383                                                if(!xIsNan<IssmDouble>(spcdata[v1]) && !xIsNan<IssmDouble>(spcdata[v2])){
     384                                                        constraints->AddObject(new SpcStatic(count+1,iomodel->numberofvertices+i+1,dof,(spcdata[v1]+spcdata[v2])/2.,analysis_type));
     385                                                        count++;
    413386                                                }
    414387                                        }

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

    r23561 r23576  
    99#include "./ModelProcessorx.h"
    1010
    11 void CreateNodes(Nodes* nodes, IoModel* iomodel,int analysis,int finite_element,int approximation){
     11void CreateNodes(Nodes* nodes, IoModel* iomodel,int analysis,int finite_element,int approximation,int* approximations){
    1212
    1313        /*Intermediaries*/
    14         int   i,j,counter,vnodes,lid=0;
    15         int   numberoffaces,elementnbv;
    16         bool *my_nodes = NULL;
    17         Node *node     = NULL;
    18 
    19         int id0 = 0;
    20         switch(finite_element){
    21                 case P1Enum:
    22                         for(i=0;i<iomodel->numberofvertices;i++){
    23                                 if(iomodel->my_vertices[i]){
    24                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    25                                 }
    26                         }
    27                         break;
    28 
    29                 case P1DGEnum:
    30                         DiscontinuousGalerkinNodesPartitioning(&my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel);
    31                         for(i=0;i<iomodel->numberofelements;i++){
    32                                 for(j=0;j<3;j++){
    33                                         if(my_nodes[3*i+j]){
    34                                                 nodes->AddObject(new Node(id0+3*i+j+1,id0+3*i+j,lid++,0,iomodel->elements[3*i+j]-1,false,iomodel,analysis,approximation));
    35 
    36                                         }
    37                                 }
    38                         }
    39                         break;
    40 
    41                 case P1bubbleEnum:
    42                         for(i=0;i<iomodel->numberofvertices;i++){
    43                                 if(iomodel->my_vertices[i]){
    44                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    45                                 }
    46                         }
    47                         for(i=0;i<iomodel->numberofelements;i++){
    48                                 if(iomodel->my_elements[i]){
    49                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation));
    50                                 }
    51                         }
    52                         break;
    53 
    54                 case P1bubblecondensedEnum:
    55                         for(i=0;i<iomodel->numberofvertices;i++){
    56                                 if(iomodel->my_vertices[i]){
    57                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    58                                 }
    59                         }
    60                         for(i=0;i<iomodel->numberofelements;i++){
    61                                 if(iomodel->my_elements[i]){
    62                                         node = new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation);
    63                                         node->HardDeactivate();
    64                                         nodes->AddObject(node);
    65                                 }
    66                         }
    67                         break;
    68 
    69                 case P1xP2Enum:
    70                         EdgesPartitioning(iomodel);
    71                         for(i=0;i<iomodel->numberofvertices;i++){
    72                                 if(iomodel->my_vertices[i]){
    73                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    74                                 }
    75                         }
    76 
    77                         counter = iomodel->numberofvertices;
    78                         for(i=0;i<iomodel->numberofedges;i++){
    79                                 if(iomodel->edges[i*3+2]==2){
    80                                         if(iomodel->my_edges[i]){
    81                                                 node = new Node(id0+iomodel->numberofvertices+i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    82                                                 nodes->AddObject(node);
    83                                         }
    84                                         counter++;
    85                                 }
    86                         }
    87                         break;
    88 
    89                 case P1xP3Enum:
    90                         EdgesPartitioning(iomodel);
    91                         for(i=0;i<iomodel->numberofvertices;i++){
    92                                 if(iomodel->my_vertices[i]){
    93                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    94                                 }
    95                         }
    96 
    97                         counter = iomodel->numberofvertices;
    98                         for(i=0;i<iomodel->numberofedges;i++){
    99                                 if(iomodel->edges[i*3+2]==2){
    100                                         if(iomodel->my_edges[i]){
    101                                                 node = new Node(id0+iomodel->numberofvertices+2*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    102                                                 nodes->AddObject(node);
    103                                                 node = new Node(id0+iomodel->numberofvertices+2*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation);
    104                                                 nodes->AddObject(node);
    105                                         }
    106                                         counter=counter+2;
    107                                 }
    108                         }
    109                         break;
    110                 case P1xP4Enum:
    111                         EdgesPartitioning(iomodel);
    112                         for(i=0;i<iomodel->numberofvertices;i++){
    113                                 if(iomodel->my_vertices[i]){
    114                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    115                                 }
    116                         }
    117                         counter = iomodel->numberofvertices;
    118                         for(i=0;i<iomodel->numberofedges;i++){
    119                                 if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
    120                                         if(iomodel->my_edges[i]){
    121                                                 node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    122                                                 nodes->AddObject(node);
    123                                                 node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation);
    124                                                 nodes->AddObject(node);
    125                                                 node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,approximation);
    126                                                 nodes->AddObject(node);
    127                                         }
    128                                         counter=counter+3;
    129                                 }
    130                         }
    131                         break;
    132 
    133                 case P2xP1Enum:
    134                         EdgesPartitioning(iomodel);
    135                         for(i=0;i<iomodel->numberofvertices;i++){
    136                                 if(iomodel->my_vertices[i]){
    137                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    138                                 }
    139                         }
    140 
    141                         counter = iomodel->numberofvertices;
    142                         for(i=0;i<iomodel->numberofedges;i++){
    143                                 if(iomodel->edges[i*3+2]!=2){
    144                                         if(iomodel->my_edges[i]){
    145                                                 node = new Node(id0+iomodel->numberofvertices+i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    146                                                 nodes->AddObject(node);
    147                                         }
    148                                         counter++;
    149                                 }
    150                         }
    151                         break;
    152 
    153                 case P2Enum:
    154                         EdgesPartitioning(iomodel);
    155                         for(i=0;i<iomodel->numberofvertices;i++){
    156                                 if(iomodel->my_vertices[i]){
    157                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    158                                 }
    159                         }
    160                         for(i=0;i<iomodel->numberofedges;i++){
    161                                 if(iomodel->my_edges[i]){
    162                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation));
    163                                 }
    164                         }
    165                         id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
    166               if(iomodel->meshelementtype==PentaEnum){
     14        int        numnodes;
     15        const int  MAXCONNECTIVITY = 5;
     16        int        element_numnodes;
     17        int        element_node_ids[40] = {0};
     18
     19        /*Get partitioning variables*/
     20        int  my_rank   = IssmComm::GetRank();
     21        int  num_procs = IssmComm::GetSize();
     22        int*     epart = iomodel->epart;
     23
     24
     25        /*Determine how many nodes we have in total*/
     26        /*Define nodes sids for each element*/
     27        if(iomodel->meshelementtype==TriaEnum){
     28                switch(finite_element){
     29                        case P1Enum:
     30                                numnodes = iomodel->numberofvertices;
     31                                break;
     32                        case P1DGEnum:
     33                                numnodes = 3*iomodel->numberofelements;
     34                                break;
     35                        case P1bubbleEnum: case P1bubblecondensedEnum:
     36                                numnodes = iomodel->numberofvertices+iomodel->numberofelements;
     37                                break;
     38                        case P2Enum:
     39                                EdgesPartitioning(iomodel);
     40                                numnodes = iomodel->numberofvertices+iomodel->numberofedges;
     41                                break;
     42                        case MINIEnum: case MINIcondensedEnum:
     43                                /*P1+ velocity (bubble statically condensed), P1 pressure*/
     44                                numnodes = (iomodel->numberofvertices+iomodel->numberofelements) + (iomodel->numberofvertices);
     45                                break;
     46                        case TaylorHoodEnum: case XTaylorHoodEnum:
     47                                /*P2 velocity, P1 pressure*/
     48                                EdgesPartitioning(iomodel);
     49                                numnodes = (iomodel->numberofvertices+iomodel->numberofedges) + (iomodel->numberofvertices);
     50                                break;
     51                        case LATaylorHoodEnum:
     52                                /*P2 velocity*/
     53                                EdgesPartitioning(iomodel);
     54                                numnodes = (iomodel->numberofvertices+iomodel->numberofedges);
     55                                break;
     56                        case CrouzeixRaviartEnum:
     57                                /*P2b velocity, P1 DG pressure*/
     58                                EdgesPartitioning(iomodel);
     59                                numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements) + (3*iomodel->numberofelements);
     60                                break;
     61                        case LACrouzeixRaviartEnum:
     62                                /*P2b velocity*/
     63                                EdgesPartitioning(iomodel);
     64                                numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements);
     65                                break;
     66                        default:
     67                                _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
     68                }
     69        }
     70        else if(iomodel->meshelementtype==PentaEnum){
     71                switch(finite_element){
     72                        case P1Enum:
     73                                numnodes = iomodel->numberofvertices;
     74                                break;
     75                        case P1bubbleEnum: case P1bubblecondensedEnum:
     76                                numnodes = iomodel->numberofvertices+iomodel->numberofelements;
     77                                break;
     78                        case P1xP2Enum:
     79                                EdgesPartitioning(iomodel);
     80                                numnodes = iomodel->numberofvertices+iomodel->numberofverticaledges;
     81                                break;
     82                        case P1xP3Enum:
     83                                EdgesPartitioning(iomodel);
     84                                numnodes = iomodel->numberofvertices+2*iomodel->numberofverticaledges;
     85                                break;
     86                        case P1xP4Enum:
     87                                EdgesPartitioning(iomodel);
     88                                numnodes = iomodel->numberofvertices+3*iomodel->numberofverticaledges;
     89                                break;
     90                        case P2xP1Enum:
     91                                EdgesPartitioning(iomodel);
     92                                numnodes = iomodel->numberofvertices+iomodel->numberofhorizontaledges;
     93                                break;
     94                        case P2Enum:
     95                                EdgesPartitioning(iomodel);
    16796                                FacesPartitioning(iomodel);
    168                                 for(i=0;i<iomodel->numberoffaces;i++){
    169                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    170                                                 if(iomodel->my_faces[i]){
    171                                                         node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,approximation);
    172                                                         nodes->AddObject(node);
     97                                numnodes = iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces;
     98                                break;
     99                        case P2xP4Enum:
     100                                EdgesPartitioning(iomodel);
     101                                FacesPartitioning(iomodel);
     102                                numnodes = iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces;
     103                                break;
     104                        case MINIEnum: case MINIcondensedEnum:
     105                                /*P1+ velocity (bubble statically condensed), P1 pressure*/
     106                                numnodes = (iomodel->numberofvertices+iomodel->numberofelements) + (iomodel->numberofvertices);
     107                                break;
     108                        case TaylorHoodEnum: case XTaylorHoodEnum:
     109                                /*P2 velocity, P1 pressure*/
     110                                EdgesPartitioning(iomodel);
     111                                FacesPartitioning(iomodel);
     112                                numnodes = (iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces) + (iomodel->numberofvertices);
     113                                break;
     114                        case OneLayerP4zEnum:
     115                                /*P2xP4 velocity, P1 pressure*/
     116                                EdgesPartitioning(iomodel);
     117                                FacesPartitioning(iomodel);
     118                                numnodes = (iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces) + (iomodel->numberofvertices);
     119                                break;
     120                        default:
     121                                _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
     122                }
     123        }
     124        else{
     125                _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
     126        }
     127
     128        /*create matrix that keeps track of all ranks that have node i, and initialize as -1 (Common to all CPUs)*/
     129        int* nodes_ranks = xNew<int>(MAXCONNECTIVITY*numnodes);
     130        for(int i=0;i<MAXCONNECTIVITY*numnodes;i++) nodes_ranks[i] = -1;
     131
     132        /*For all nodes, count how many cpus have node i (initialize with 0)*/
     133        int* nodes_proc_count = xNewZeroInit<int>(numnodes);
     134
     135        /*Create vector of size total numnodes, initialized with -1, that will keep track of local ids*/
     136        int* nodes_lids  = xNew<int>(numnodes);
     137        for(int i=0;i<numnodes;i++) nodes_lids[i] = -1;
     138
     139        /*Create vector of approximation per node (used for FS: vel or pressure)*/
     140        int* nodes_approx = xNew<int>(numnodes);
     141        if(approximations){
     142                for(int i=0;i<numnodes;i++) nodes_approx[i] = approximations[i];
     143        }
     144        else{
     145                for(int i=0;i<numnodes;i++) nodes_approx[i] = approximation;
     146        }
     147
     148        /*Go through all elements and mark all vertices for all partitions*/
     149        int  lid = 0;
     150        for(int i=0;i<iomodel->numberofelements;i++){
     151
     152                /*Define nodes sids for each element*/
     153                if(iomodel->meshelementtype==TriaEnum){
     154                        switch(finite_element){
     155                                case P1Enum:
     156                                        element_numnodes=3;
     157                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     158                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     159                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     160                                        break;
     161                                case P1bubbleEnum: case P1bubblecondensedEnum:
     162                                        element_numnodes=4;
     163                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     164                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     165                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     166                                        element_node_ids[3]=iomodel->numberofvertices+i;
     167                                        break;
     168                                case P1DGEnum:
     169                                        element_numnodes=3;
     170                                        element_node_ids[0]=3*i+0;
     171                                        element_node_ids[1]=3*i+1;
     172                                        element_node_ids[2]=3*i+2;
     173                                        break;
     174                                case P2Enum:
     175                                        element_numnodes = 6;
     176                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     177                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     178                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     179                                        element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
     180                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
     181                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
     182                                        break;
     183                                case MINIEnum: case MINIcondensedEnum:
     184                                        element_numnodes = 4+3;
     185                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     186                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     187                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     188                                        element_node_ids[3]=iomodel->numberofvertices+i;
     189                                        for(int n=0;n<4;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     190                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+0]-1;
     191                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+1]-1;
     192                                        element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*i+2]-1;
     193                                        for(int n=4;n<7;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
     194                                        break;
     195                                case TaylorHoodEnum: case XTaylorHoodEnum:
     196                                        element_numnodes = 6+3;
     197                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     198                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     199                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     200                                        element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
     201                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
     202                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
     203                                        for(int n=0;n<6;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     204                                        element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+0]-1;
     205                                        element_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+1]-1;
     206                                        element_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*i+2]-1;
     207                                        for(int n=6;n<9;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
     208                                        break;
     209                                case LATaylorHoodEnum:
     210                                        element_numnodes = 6;
     211                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     212                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     213                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     214                                        element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
     215                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
     216                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
     217                                        for(int n=0;n<6;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     218                                        break;
     219                                case CrouzeixRaviartEnum:
     220                                        element_numnodes = 7+3;
     221                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     222                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     223                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     224                                        element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
     225                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
     226                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
     227                                        element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+i;
     228                                        for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     229                                        element_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+0;
     230                                        element_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+1;
     231                                        element_node_ids[9]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*i+2;
     232                                        for(int n=7;n<10;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
     233                                        break;
     234                                case LACrouzeixRaviartEnum:
     235                                        element_numnodes = 7;
     236                                        element_node_ids[0]=iomodel->elements[3*i+0]-1;
     237                                        element_node_ids[1]=iomodel->elements[3*i+1]-1;
     238                                        element_node_ids[2]=iomodel->elements[3*i+2]-1;
     239                                        element_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+0];
     240                                        element_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+1];
     241                                        element_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*i+2];
     242                                        element_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+i;
     243                                        for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     244                                        break;
     245                                default:
     246                                        _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
     247                        }
     248                }
     249                else if(iomodel->meshelementtype==PentaEnum){
     250                        switch(finite_element){
     251                                case P1Enum:
     252                                        element_numnodes=6;
     253                                        element_node_ids[0]=iomodel->elements[6*i+0]-1;
     254                                        element_node_ids[1]=iomodel->elements[6*i+1]-1;
     255                                        element_node_ids[2]=iomodel->elements[6*i+2]-1;
     256                                        element_node_ids[3]=iomodel->elements[6*i+3]-1;
     257                                        element_node_ids[4]=iomodel->elements[6*i+4]-1;
     258                                        element_node_ids[5]=iomodel->elements[6*i+5]-1;
     259                                        break;
     260                                case P1bubbleEnum: case P1bubblecondensedEnum:
     261                                        element_numnodes=7;
     262                                        element_node_ids[0]=iomodel->elements[6*i+0]-1;
     263                                        element_node_ids[1]=iomodel->elements[6*i+1]-1;
     264                                        element_node_ids[2]=iomodel->elements[6*i+2]-1;
     265                                        element_node_ids[3]=iomodel->elements[6*i+3]-1;
     266                                        element_node_ids[4]=iomodel->elements[6*i+4]-1;
     267                                        element_node_ids[5]=iomodel->elements[6*i+5]-1;
     268                                        element_node_ids[6]=iomodel->numberofvertices+i;
     269                                        break;
     270                                case P1xP2Enum:
     271                                        element_numnodes=9;
     272                                        element_node_ids[0]=iomodel->elements[6*i+0]-1;
     273                                        element_node_ids[1]=iomodel->elements[6*i+1]-1;
     274                                        element_node_ids[2]=iomodel->elements[6*i+2]-1;
     275                                        element_node_ids[3]=iomodel->elements[6*i+3]-1;
     276                                        element_node_ids[4]=iomodel->elements[6*i+4]-1;
     277                                        element_node_ids[5]=iomodel->elements[6*i+5]-1;
     278                                        element_node_ids[6]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+0];
     279                                        element_node_ids[7]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+1];
     280                                        element_node_ids[8]=iomodel->numberofvertices+iomodel->elementtoverticaledgeconnectivity[3*i+2];
     281                                        break;
     282                                case P1xP3Enum:
     283                                        element_numnodes=12;
     284                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
     285                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
     286                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
     287                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
     288                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
     289                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
     290                                        element_node_ids[ 6]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+0;
     291                                        element_node_ids[ 7]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+0;
     292                                        element_node_ids[ 8]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+0;
     293                                        element_node_ids[ 9]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1;
     294                                        element_node_ids[10]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1;
     295                                        element_node_ids[11]=iomodel->numberofvertices+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1;
     296                                        break;
     297                                case P1xP4Enum:
     298                                        element_numnodes=15;
     299                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/
     300                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/
     301                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/
     302                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/
     303                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/
     304                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/
     305                                        element_node_ids[ 6]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+0; /*mid vertical edge 1*/
     306                                        element_node_ids[ 7]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+0; /*mid vertical edge 2*/
     307                                        element_node_ids[ 8]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+0; /*mid vertical edge 3*/
     308                                        element_node_ids[ 9]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 1/4 vertical edge 1*/
     309                                        element_node_ids[10]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 1/4 vertical edge 2*/
     310                                        element_node_ids[11]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 1/4 vertical edge 3*/
     311                                        element_node_ids[12]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+0]+2; /* 3/4 vertical edge 1*/
     312                                        element_node_ids[13]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+1]+2; /* 3/4 vertical edge 2*/
     313                                        element_node_ids[14]=iomodel->numberofvertices+3*iomodel->elementtoverticaledgeconnectivity[3*i+2]+2; /* 3/4 vertical edge 3*/
     314                                        break;
     315                                case P2xP1Enum:
     316                                        element_numnodes=12;
     317                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
     318                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
     319                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
     320                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
     321                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
     322                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
     323                                        element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+0];
     324                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+1];
     325                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+2];
     326                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+3];
     327                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+4];
     328                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtohorizontaledgeconnectivity[6*i+5];
     329                                        break;
     330                                case P2Enum:
     331                                        element_numnodes = 18;
     332                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
     333                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
     334                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
     335                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
     336                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
     337                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
     338                                        element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0];
     339                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1];
     340                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2];
     341                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3];
     342                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4];
     343                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5];
     344                                        element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6];
     345                                        element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7];
     346                                        element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8];
     347                                        element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+0];
     348                                        element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+1];
     349                                        element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+2];
     350                                        break;
     351                                case P2xP4Enum:
     352                                        element_numnodes = 30;
     353                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/
     354                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/
     355                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/
     356                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/
     357                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/
     358                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/
     359                                        element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; /*mid vertical edge 1*/
     360                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; /*mid vertical edge 2*/
     361                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; /*mid vertical edge 3*/
     362                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; /*mid basal edge 1*/
     363                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; /*mid basal edge 2*/
     364                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; /*mid basal edge 3*/
     365                                        element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; /*mid top edge 1*/
     366                                        element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; /*mid top edge 2*/
     367                                        element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; /*mid top edge 3*/
     368                                        element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]; /* 1/4 vertical edge 1*/
     369                                        element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]; /* 1/4 vertical edge 2*/
     370                                        element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]; /* 1/4 vertical edge 3*/
     371                                        element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 3/4 vertical edge 1*/
     372                                        element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 3/4 vertical edge 2*/
     373                                        element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 3/4 vertical edge 3*/
     374                                        element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+0; /* 1/4 vertical face 1*/
     375                                        element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+0; /* 1/4 vertical face 2*/
     376                                        element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+0; /* 1/4 vertical face 3*/
     377                                        element_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+1; /* 2/4 vertical face 1*/
     378                                        element_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+1; /* 2/4 vertical face 2*/
     379                                        element_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+1; /* 2/4 vertical face 3*/
     380                                        element_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+2; /* 3/4 vertical face 1*/
     381                                        element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/
     382                                        element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/
     383                                        break;
     384                                case MINIEnum: case MINIcondensedEnum:
     385                                        element_numnodes = 7+6;
     386                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
     387                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
     388                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
     389                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
     390                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
     391                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
     392                                        element_node_ids[ 6]=iomodel->numberofvertices+i;
     393                                        if(!approximations) for(int n=0;n<7;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     394                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+0]-1;
     395                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+1]-1;
     396                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+2]-1;
     397                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+3]-1;
     398                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+4]-1;
     399                                        element_node_ids[12]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[6*i+5]-1;
     400                                        if(!approximations) for(int n=7;n<13;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
     401                                        break;
     402                                case TaylorHoodEnum: case XTaylorHoodEnum:
     403                                        element_numnodes = 18+6;
     404                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1;
     405                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1;
     406                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1;
     407                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1;
     408                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1;
     409                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1;
     410                                        element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0];
     411                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1];
     412                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2];
     413                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3];
     414                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4];
     415                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5];
     416                                        element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6];
     417                                        element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7];
     418                                        element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8];
     419                                        element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+0];
     420                                        element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+1];
     421                                        element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtoverticalfaceconnectivity[3*i+2];
     422                                        for(int n=0;n<18;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     423                                        element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+0]-1;
     424                                        element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+1]-1;
     425                                        element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+2]-1;
     426                                        element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+3]-1;
     427                                        element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+4]-1;
     428                                        element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofverticalfaces+iomodel->elements[6*i+5]-1;
     429                                        for(int n=18;n<24;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
     430                                        break;
     431                                case OneLayerP4zEnum:
     432                                        element_numnodes = 30+6;
     433                                        element_node_ids[ 0]=iomodel->elements[6*i+0]-1; /*Vertex 1*/
     434                                        element_node_ids[ 1]=iomodel->elements[6*i+1]-1; /*Vertex 2*/
     435                                        element_node_ids[ 2]=iomodel->elements[6*i+2]-1; /*Vertex 3*/
     436                                        element_node_ids[ 3]=iomodel->elements[6*i+3]-1; /*Vertex 4*/
     437                                        element_node_ids[ 4]=iomodel->elements[6*i+4]-1; /*Vertex 5*/
     438                                        element_node_ids[ 5]=iomodel->elements[6*i+5]-1; /*Vertex 6*/
     439                                        element_node_ids[ 6]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+0]; /*mid vertical edge 1*/
     440                                        element_node_ids[ 7]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+1]; /*mid vertical edge 2*/
     441                                        element_node_ids[ 8]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+2]; /*mid vertical edge 3*/
     442                                        element_node_ids[ 9]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+3]; /*mid basal edge 1*/
     443                                        element_node_ids[10]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+4]; /*mid basal edge 2*/
     444                                        element_node_ids[11]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+5]; /*mid basal edge 3*/
     445                                        element_node_ids[12]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+6]; /*mid top edge 1*/
     446                                        element_node_ids[13]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+7]; /*mid top edge 2*/
     447                                        element_node_ids[14]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*i+8]; /*mid top edge 3*/
     448                                        element_node_ids[15]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]; /* 1/4 vertical edge 1*/
     449                                        element_node_ids[16]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]; /* 1/4 vertical edge 2*/
     450                                        element_node_ids[17]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]; /* 1/4 vertical edge 3*/
     451                                        element_node_ids[18]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+0]+1; /* 3/4 vertical edge 1*/
     452                                        element_node_ids[19]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+1]+1; /* 3/4 vertical edge 2*/
     453                                        element_node_ids[20]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->elementtoverticaledgeconnectivity[3*i+2]+1; /* 3/4 vertical edge 3*/
     454                                        element_node_ids[21]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+0; /* 1/4 vertical face 1*/
     455                                        element_node_ids[22]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+0; /* 1/4 vertical face 2*/
     456                                        element_node_ids[23]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+0; /* 1/4 vertical face 3*/
     457                                        element_node_ids[24]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+1; /* 2/4 vertical face 1*/
     458                                        element_node_ids[25]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+1; /* 2/4 vertical face 2*/
     459                                        element_node_ids[26]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+1; /* 2/4 vertical face 3*/
     460                                        element_node_ids[27]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+0]+2; /* 3/4 vertical face 1*/
     461                                        element_node_ids[28]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+1]+2; /* 3/4 vertical face 2*/
     462                                        element_node_ids[29]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->elementtoverticalfaceconnectivity[3*i+2]+2; /* 3/4 vertical face 3*/
     463                                        for(int n=0;n<30;n++) nodes_approx[element_node_ids[n]] = FSvelocityEnum;
     464                                        element_node_ids[30]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+0]-1;
     465                                        element_node_ids[31]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+1]-1;
     466                                        element_node_ids[32]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+2]-1;
     467                                        element_node_ids[33]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+3]-1;
     468                                        element_node_ids[34]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+4]-1;
     469                                        element_node_ids[35]=iomodel->numberofvertices+iomodel->numberofedges+2*iomodel->numberofverticaledges+3*iomodel->numberofverticalfaces+iomodel->elements[6*i+5]-1;
     470                                        for(int n=30;n<36;n++) nodes_approx[element_node_ids[n]] = FSpressureEnum;
     471                                        break;
     472                                default:
     473                                        _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
     474                        }
     475                }
     476                else{
     477                        _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
     478                }
     479
     480                for(int j=0;j<element_numnodes;j++){
     481                        int nid = element_node_ids[j]; _assert_(nid<numnodes);
     482
     483                        /*See if it has already been marked*/
     484                        bool found = false;
     485                        for(int k=0;k<nodes_proc_count[nid];k++){
     486                                if(nodes_ranks[MAXCONNECTIVITY*nid+k] == epart[i]){
     487                                        found = true;
     488                                        break;
     489                                }
     490                        }
     491
     492                        /*On go below if this vertex has not been seen yet in this partition*/
     493                        if(!found){
     494                                /*This rank has not been marked for this vertex just yet so go ahead and mark it*/
     495                                nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = epart[i];
     496                                nodes_proc_count[nid]++;
     497
     498                                /*Keep track of local ids!*/
     499                                if(epart[i]==my_rank){
     500                                        nodes_lids[nid] = lid;
     501                                        lid++;
     502                                }
     503
     504                                /*Make sure we don't go too far in the table*/
     505                                if(nodes_proc_count[nid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
     506                        }
     507                }
     508        }
     509        if(finite_element==P1DGEnum){/* Special case for DG...{{{*/
     510                int node_list[4];
     511                if(!iomodel->domaintype==Domain2DhorizontalEnum) _error_("not implemented yet");
     512                CreateEdges(iomodel);
     513                CreateFaces(iomodel);
     514                for(int i=0;i<iomodel->numberoffaces;i++){
     515                        int e1=iomodel->faces[4*i+2]-1; //faces are [node1 node2 elem1 elem2]
     516                        int e2=iomodel->faces[4*i+3]-1; //faces are [node1 node2 elem1 elem2]
     517                        if(e2!=-2){
     518                                if(epart[e1]!=epart[e2]){
     519                                        int i1=iomodel->faces[4*i+0];
     520                                        int i2=iomodel->faces[4*i+1];
     521                                        int pos=-1;
     522                                        for(int j=0;j<3;j++) if(iomodel->elements[3*e2+j]==i1) pos=j;
     523                                        if(     pos==0){ node_list[0] = e2*3+0; node_list[1] = e2*3+2;}
     524                                        else if(pos==1){ node_list[0] = e2*3+1; node_list[1] = e2*3+0;}
     525                                        else if(pos==2){ node_list[0] = e2*3+2; node_list[1] = e2*3+1;}
     526                                        else _error_("not supposed to happen");
     527                                        pos=-1;
     528                                        for(int j=0;j<3;j++) if(iomodel->elements[3*e1+j]==i1) pos=j;
     529                                        if(     pos==0){ node_list[2] = e1*3+0; node_list[3] = e1*3+1;}
     530                                        else if(pos==1){ node_list[2] = e1*3+1; node_list[3] = e1*3+2;}
     531                                        else if(pos==2){ node_list[2] = e1*3+2; node_list[3] = e1*3+0;}
     532                                        else _error_("not supposed to happen");
     533                                        for(int j=0;j<4;j++){
     534                                                int  nid = node_list[j];
     535                                                bool found = false;
     536                                                for(int k=0;k<nodes_proc_count[nid];k++){
     537                                                        if(nodes_ranks[MAXCONNECTIVITY*nid+k] == epart[e1]){
     538                                                                found = true;
     539                                                                break;
     540                                                        }
     541                                                }
     542                                                if(!found){
     543                                                        nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = epart[e1];
     544                                                        nodes_proc_count[nid]++;
     545                                                        if(epart[e1]==my_rank){
     546                                                                nodes_lids[nid] = lid;
     547                                                                lid++;
     548                                                        }
     549                                                        if(nodes_proc_count[nid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
    173550                                                }
    174551                                        }
    175                                         else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    176                                                 /*Nothing*/
    177                                         }
    178                                         else{
    179                                                 _error_("not supported");
    180                                         }
    181                                 }
    182                         }
    183                         break;
    184                 case P2bubbleEnum:
    185                         EdgesPartitioning(iomodel);
    186                         for(i=0;i<iomodel->numberofvertices;i++){
    187                                 if(iomodel->my_vertices[i]){
    188                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    189                                 }
    190                         }
    191                         for(i=0;i<iomodel->numberofedges;i++){
    192                                 if(iomodel->my_edges[i]){
    193                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,approximation));
    194                                 }
    195                         }
    196                         id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
    197                         if(iomodel->meshelementtype==PentaEnum){
    198                                 FacesPartitioning(iomodel);
    199                                 for(i=0;i<iomodel->numberoffaces;i++){
    200                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    201                                                 if(iomodel->my_faces[i]){
    202                                                         node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,approximation);
    203                                                         nodes->AddObject(node);
     552                                        for(int j=0;j<4;j++){
     553                                                int  nid = node_list[j];
     554                                                bool found = false;
     555                                                for(int k=0;k<nodes_proc_count[nid];k++){
     556                                                        if(nodes_ranks[MAXCONNECTIVITY*nid+k] == epart[e2]){
     557                                                                found = true;
     558                                                                break;
     559                                                        }
     560                                                }
     561                                                if(!found){
     562                                                        nodes_ranks[MAXCONNECTIVITY*nid+nodes_proc_count[nid]] = epart[e2];
     563                                                        nodes_proc_count[nid]++;
     564                                                        if(epart[e2]==my_rank){
     565                                                                nodes_lids[nid] = lid;
     566                                                                lid++;
     567                                                        }
     568                                                        if(nodes_proc_count[nid]>MAXCONNECTIVITY) _error_("need to increase MAXCONNECTIVITY (this is hard coded, contact ISSM developer)");
    204569                                                }
    205570                                        }
    206                                         else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    207                                                 /*Nothing*/
     571                                }
     572                        }
     573                }
     574        }/*}}}*/
     575        //if(my_rank==0) printarray(nodes_ranks,numnodes,MAXCONNECTIVITY);
     576
     577        /*Now, Count how many clones we have with other partitions*/
     578        int*  common_send = xNew<int>(num_procs);
     579        int*  common_recv = xNew<int>(num_procs);
     580        int** common_send_ids = xNew<int*>(num_procs);
     581        int** common_recv_ids = xNew<int*>(num_procs);
     582
     583        /*First step: allocate, Step 2: populate*/
     584        for(int step=0;step<2;step++){
     585                if(step==1){
     586                        /*Allocate send and receive arrays of ids now*/
     587                        for(int i=0;i<num_procs;i++){
     588                                _assert_(common_send[i]>=0 && common_recv[i]>=0);
     589                                common_send_ids[i] = xNew<int>(common_send[i]);
     590                                common_recv_ids[i] = xNew<int>(common_recv[i]);
     591                        }
     592                }
     593                /*Re/Initialize counters to 0*/
     594                for(int i=0;i<num_procs;i++){
     595                        common_recv[i]=0;
     596                        common_send[i]=0;
     597                }
     598                /*Go through table and find clones/masters etc*/
     599                for(int i=0;i<numnodes;i++){
     600                        /*If we did not find this vertex in our current partition, go to next vertex*/
     601                        if(nodes_lids[i] == -1) continue;
     602                        /*Find in what column this rank belongs*/
     603                        int col = -1;
     604                        for(int j=0;j<MAXCONNECTIVITY;j++){
     605                                if(nodes_ranks[MAXCONNECTIVITY*i+j] == my_rank){
     606                                        col = j;
     607                                        break;
     608                                }
     609                        }
     610                        _assert_(col!=-1);
     611
     612                        /*If col==0, it is either not on boundary, or a master*/
     613                        if(col==0){
     614                                /*1. is this vertex on the boundary? Skip if not*/
     615                                if(nodes_ranks[MAXCONNECTIVITY*i+col+1]==-1){
     616                                        continue;
     617                                }
     618                                else{
     619                                        for(int j=1;j<nodes_proc_count[i];j++){
     620                                                _assert_(nodes_ranks[MAXCONNECTIVITY*i+j]>=0);
     621                                                int rank = nodes_ranks[MAXCONNECTIVITY*i+j];
     622                                                if(step==1){
     623                                                        common_send_ids[rank][common_send[rank]] = nodes_lids[i];
     624                                                }
     625                                                common_send[rank]++;
    208626                                        }
    209                                         else{
    210                                                 _error_("not supported");
    211                                         }
    212                                 }
    213                                 id0 = id0+iomodel->numberoffaces;
    214                         }
    215                         for(i=0;i<iomodel->numberofelements;i++){
    216                                 if(iomodel->my_elements[i]){
    217                                         nodes->AddObject(new Node(id0+i+1,id0+i,lid++,0,0,false,iomodel,analysis,approximation));
    218                                 }
    219                         }
    220                         break;
    221                 case P2xP4Enum:
    222                         EdgesPartitioning(iomodel);
    223                         FacesPartitioning(iomodel);
    224                         for(i=0;i<iomodel->numberofvertices;i++){
    225                                 if(iomodel->my_vertices[i]){
    226                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,approximation));
    227                                 }
    228                         }
    229                         counter = iomodel->numberofvertices;
    230                         for(i=0;i<iomodel->numberofedges;i++){
    231                                 if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
    232                                         if(iomodel->my_edges[i]){
    233                                                 node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    234                                                 nodes->AddObject(node);
    235                                                 node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation);
    236                                                 nodes->AddObject(node);
    237                                                 node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,approximation);
    238                                                 nodes->AddObject(node);
    239                                         }
    240                                         counter=counter+3;
    241                                 }
    242                                 else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
    243                                         if(iomodel->my_edges[i]){
    244                                                 node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    245                                                 nodes->AddObject(node);
    246                                         }
    247                                         counter=counter+1;
    248                                 }
    249                                 else{
    250                                         _error_("not supported");
    251                                 }
    252                         }
    253                         id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges;
    254                         for(i=0;i<iomodel->numberoffaces;i++){
    255                                 if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    256                                         if(iomodel->my_faces[i]){
    257                                                 node = new Node(id0+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,approximation);
    258                                                 nodes->AddObject(node);
    259                                                 node = new Node(id0+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,approximation);
    260                                                 nodes->AddObject(node);
    261                                                 node = new Node(id0+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,approximation);
    262                                                 nodes->AddObject(node);
    263                                         }
    264                                         counter=counter+3;
    265                                 }
    266                                 else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    267                                         /*Nothing*/
    268                                 }
    269                                 else{
    270                                         _error_("not supported");
    271                                 }
    272                         }
    273                         break;
    274 
    275                 /*Stokes elements*/
    276                 case P1P1Enum:
    277                         _assert_(approximation==FSApproximationEnum);
    278                         /*P1 velocity*/
    279                         for(i=0;i<iomodel->numberofvertices;i++){
    280                                 if(iomodel->my_vertices[i]){
    281                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    282                                 }
    283                         }
    284                         /*P1 pressure*/
    285                         vnodes = id0+iomodel->numberofvertices;
    286                         for(i=0;i<iomodel->numberofvertices;i++){
    287                                 if(iomodel->my_vertices[i]){
    288                                         nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum));
    289                                 }
    290                         }
    291                         break;
    292                 case P1P1GLSEnum:
    293                         _assert_(approximation==FSApproximationEnum);
    294                         /*P1 velocity*/
    295                         for(i=0;i<iomodel->numberofvertices;i++){
    296                                 if(iomodel->my_vertices[i]){
    297                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    298                                 }
    299                         }
    300                         /*P1 pressure*/
    301                         vnodes = id0+iomodel->numberofvertices;
    302                         for(i=0;i<iomodel->numberofvertices;i++){
    303                                 if(iomodel->my_vertices[i]){
    304                                         nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum));
    305                                 }
    306                         }
    307                         break;
    308                 case MINIcondensedEnum:
    309                         _assert_(approximation==FSApproximationEnum);
    310                         /*P1+ velocity (bubble statically condensed)*/
    311                         for(i=0;i<iomodel->numberofvertices;i++){
    312                                 if(iomodel->my_vertices[i]){
    313                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    314                                 }
    315                         }
    316                         for(i=0;i<iomodel->numberofelements;i++){
    317                                 if(iomodel->my_elements[i]){
    318                                         node = new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
     627                                }
     628                        }
     629                        else{
     630                                /*3. It is a slave, record that we need to receive for this cpu*/
     631                                int rank = nodes_ranks[MAXCONNECTIVITY*i+0];
     632                                if(step==1){
     633                                        common_recv_ids[rank][common_recv[rank]] = nodes_lids[i];
     634                                }
     635                                common_recv[rank]++;
     636                        }
     637                }
     638        }
     639        xDelete<int>(nodes_proc_count);
     640
     641        /*Final step: prepare pids (parallel ids), first count number of masters for each proc*/
     642        int*  num_masters = xNewZeroInit<int>(num_procs);
     643        for(int i=0;i<numnodes;i++){
     644                num_masters[nodes_ranks[MAXCONNECTIVITY*i+0]]++;
     645        }
     646        /*Now count offsets for each procs (by taking the sum of masters of procs<my_rank)*/
     647        int* rank_offsets = xNewZeroInit<int>(num_procs);
     648        for(int i=0;i<num_procs;i++){
     649                for(int j=i+1;j<num_procs;j++) rank_offsets[i]+=num_masters[j];
     650        }
     651        xDelete<int>(num_masters);
     652        /*Now create pids*/
     653        int* nodes_pids = xNew<int>(numnodes);
     654        int* rank_counters = xNewZeroInit<int>(numnodes);
     655        for(int i=0;i<numnodes;i++){
     656                int rank_master = nodes_ranks[MAXCONNECTIVITY*i+0];
     657                nodes_pids[i] = rank_offsets[rank_master]+rank_counters[rank_master];
     658                rank_counters[rank_master]++;
     659        }
     660        xDelete<int>(rank_counters);
     661        xDelete<int>(rank_offsets);
     662
     663        /*Go ahead and create vertices now that we have all we need*/
     664        for(int i=0;i<numnodes;i++){
     665                if(nodes_lids[i]!=-1){
     666                        bool isclone = (nodes_ranks[MAXCONNECTIVITY*i+0]!=my_rank);
     667                        int io_index = 0;
     668                        if(i<iomodel->numberofvertices) io_index = i;
     669                        Node* node=new Node(i+1,i,nodes_lids[i],nodes_pids[i],io_index,isclone,iomodel,analysis,nodes_approx[i]);
     670                        if(finite_element==MINIcondensedEnum || finite_element==P1bubblecondensedEnum){
     671                                /*Bubble function is collapsed, needs to constrain it, maybe this is not the best place to do this, but that's life!*/
     672                                if(i>=iomodel->numberofvertices && i<iomodel->numberofvertices+iomodel->numberofelements){
    319673                                        node->HardDeactivate();
    320                                         nodes->AddObject(node);
    321                                 }
    322                         }
    323                         /*P1 pressure*/
    324                         vnodes = id0+iomodel->numberofvertices+iomodel->numberofelements;
    325                         for(i=0;i<iomodel->numberofvertices;i++){
    326                                 if(iomodel->my_vertices[i]){
    327                                         nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum));
    328                                 }
    329                         }
    330                         break;
    331                 case MINIEnum:
    332                         _assert_(approximation==FSApproximationEnum);
    333                         /*P1+ velocity*/
    334                         for(i=0;i<iomodel->numberofvertices;i++){
    335                                 if(iomodel->my_vertices[i]){
    336                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    337                                 }
    338                         }
    339                         for(i=0;i<iomodel->numberofelements;i++){
    340                                 if(iomodel->my_elements[i]){
    341                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    342                                 }
    343                         }
    344                         /*P1 pressure*/
    345                         vnodes = id0+iomodel->numberofvertices+iomodel->numberofelements;
    346                         for(i=0;i<iomodel->numberofvertices;i++){
    347                                 if(iomodel->my_vertices[i]){
    348                                         nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum));
    349                                 }
    350                         }
    351                         break;
    352                 case TaylorHoodEnum:
    353                 case XTaylorHoodEnum:
    354                         _assert_(approximation==FSApproximationEnum);
    355                         /*P2 velocity*/
    356                         EdgesPartitioning(iomodel);
    357                         for(i=0;i<iomodel->numberofvertices;i++){
    358                                 if(iomodel->my_vertices[i]){
    359                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    360                                 }
    361                         }
    362                         for(i=0;i<iomodel->numberofedges;i++){
    363                                 if(iomodel->my_edges[i]){
    364                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    365                                 }
    366                         }
    367                         id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
    368               if(iomodel->meshelementtype==PentaEnum){
    369                                 FacesPartitioning(iomodel);
    370                                 for(i=0;i<iomodel->numberoffaces;i++){
    371                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    372                                                 if(iomodel->my_faces[i]){
    373                                                         node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    374                                                         nodes->AddObject(node);
    375                                                 }
    376                                         }
    377                                         else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    378                                                 /*Nothing*/
    379                                         }
    380                                         else{
    381                                                 _error_("not supported");
    382                                         }
    383                                 }
    384                         }
    385 
    386                         /*P1 pressure*/
    387               if(iomodel->meshelementtype==PentaEnum){
    388                                 numberoffaces=iomodel->numberoffaces;
    389                         }
    390                         else{
    391                                 numberoffaces=0;
    392                         }
    393                         vnodes = id0+numberoffaces;
    394                         for(i=0;i<iomodel->numberofvertices;i++){
    395                                 if(iomodel->my_vertices[i]){
    396                                         nodes->AddObject(new Node(vnodes+i+1,iomodel->numberofvertices+iomodel->numberofedges+numberoffaces+i,lid++,0,i,false,iomodel,analysis,FSpressureEnum));
    397                                 }
    398                         }
    399                         break;
    400                 case LATaylorHoodEnum:
    401                         _assert_(approximation==FSApproximationEnum);
    402                         /*P2 velocity*/
    403                         EdgesPartitioning(iomodel);
    404                         for(i=0;i<iomodel->numberofvertices;i++){
    405                                 if(iomodel->my_vertices[i]){
    406                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    407                                 }
    408                         }
    409                         for(i=0;i<iomodel->numberofedges;i++){
    410                                 if(iomodel->my_edges[i]){
    411                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    412                                 }
    413                         }
    414                         id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
    415               if(iomodel->meshelementtype==PentaEnum){
    416                                 FacesPartitioning(iomodel);
    417                                 for(i=0;i<iomodel->numberoffaces;i++){
    418                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    419                                                 if(iomodel->my_faces[i]){
    420                                                         node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    421                                                         nodes->AddObject(node);
    422                                                 }
    423                                         }
    424                                         else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    425                                                 /*Nothing*/
    426                                         }
    427                                         else{
    428                                                 _error_("not supported");
    429                                         }
    430                                 }
    431                         }
    432 
    433                         /*No pressure*/
    434                         break;
    435                 case OneLayerP4zEnum:
    436                         _assert_(approximation==FSApproximationEnum);
    437                         EdgesPartitioning(iomodel);
    438                         FacesPartitioning(iomodel);
    439                         /*P2xP4 velocity*/
    440                         for(i=0;i<iomodel->numberofvertices;i++){
    441                                 if(iomodel->my_vertices[i]){
    442                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    443                                 }
    444                         }
    445                         counter = iomodel->numberofvertices;
    446                         for(i=0;i<iomodel->numberofedges;i++){
    447                                 if(iomodel->edges[i*3+2]==2){/*Vertical edges*/
    448                                         if(iomodel->my_edges[i]){
    449                                                 node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    450                                                 nodes->AddObject(node);
    451                                                 node = new Node(id0+iomodel->numberofvertices+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    452                                                 nodes->AddObject(node);
    453                                                 node = new Node(id0+iomodel->numberofvertices+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    454                                                 nodes->AddObject(node);
    455                                         }
    456                                         counter=counter+3;
    457                                 }
    458                                 else if(iomodel->edges[i*3+2]==1){/*Horizontal edges*/
    459                                         if(iomodel->my_edges[i]){
    460                                                 node = new Node(id0+iomodel->numberofvertices+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    461                                                 nodes->AddObject(node);
    462                                         }
    463                                         counter=counter+1;
    464                                 }
    465                                 else{
    466                                         _error_("not supported");
    467                                 }
    468                         }
    469                         id0 = id0+iomodel->numberofvertices+3*iomodel->numberofedges;
    470                         for(i=0;i<iomodel->numberoffaces;i++){
    471                                 if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    472                                         if(iomodel->my_faces[i]){
    473                                                 node = new Node(id0+3*i+1,counter+1,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    474                                                 nodes->AddObject(node);
    475                                                 node = new Node(id0+3*i+2,counter+2,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    476                                                 nodes->AddObject(node);
    477                                                 node = new Node(id0+3*i+3,counter+3,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    478                                                 nodes->AddObject(node);
    479                                         }
    480                                         counter=counter+3;
    481                                 }
    482                                 else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    483                                         /*Nothing*/
    484                                 }
    485                                 else{
    486                                         _error_("not supported");
    487                                 }
    488                         }
    489 
    490                         /*P1 pressure*/
    491                         vnodes = id0+3*iomodel->numberoffaces;
    492                         for(i=0;i<iomodel->numberofvertices;i++){
    493                                 if(iomodel->my_vertices[i]){
    494                                         nodes->AddObject(new Node(vnodes+i+1,counter+1,lid++,0,i,false,iomodel,analysis,FSpressureEnum));
    495                                 }
    496                                 counter++;
    497                         }
    498                         break;
    499                 case CrouzeixRaviartEnum:
    500                         _assert_(approximation==FSApproximationEnum);
    501                         /*P2b velocity*/
    502                         EdgesPartitioning(iomodel);
    503                         for(i=0;i<iomodel->numberofvertices;i++){
    504                                 if(iomodel->my_vertices[i]){
    505                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    506                                 }
    507                         }
    508                         for(i=0;i<iomodel->numberofedges;i++){
    509                                 if(iomodel->my_edges[i]){
    510                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    511                                 }
    512                         }
    513                         id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
    514                         if(iomodel->meshelementtype==PentaEnum){
    515                                 FacesPartitioning(iomodel);
    516                                 for(i=0;i<iomodel->numberoffaces;i++){
    517                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    518                                                 if(iomodel->my_faces[i]){
    519                                                         node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    520                                                         nodes->AddObject(node);
    521                                                 }
    522                                         }
    523                                         else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    524                                                 /*Nothing*/
    525                                         }
    526                                         else{
    527                                                 _error_("not supported");
    528                                         }
    529                                 }
    530                                 id0 = id0+iomodel->numberoffaces;
    531                         }
    532                         for(i=0;i<iomodel->numberofelements;i++){
    533                                 if(iomodel->my_elements[i]){
    534                                         nodes->AddObject(new Node(id0+i+1,id0+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    535                                 }
    536                         }
    537 
    538                         /*P1 DG pressure*/
    539                         vnodes = id0+iomodel->numberofelements;
    540                         switch(iomodel->meshelementtype){
    541                                 case TriaEnum:  elementnbv = 3; break;
    542                                 case TetraEnum: elementnbv = 4; break;
    543                                 case PentaEnum: elementnbv = 6; break;
    544                                 default:        _error_("mesh dimension not supported yet");
    545                         }
    546                         for(i=0;i<iomodel->numberofelements;i++){
    547                                 if(iomodel->my_elements[i]){
    548                                         for(j=0;j<elementnbv;j++){
    549                                                 nodes->AddObject(new Node(vnodes+elementnbv*i+j+1,vnodes+elementnbv*i+j,lid++,0,iomodel->elements[+elementnbv*i+j]-1,false,iomodel,analysis,FSpressureEnum));
    550 
    551                                         }
    552                                 }
    553                         }
    554                         break;
    555                 case LACrouzeixRaviartEnum:
    556                         _assert_(approximation==FSApproximationEnum);
    557                         /*P2b velocity*/
    558                         EdgesPartitioning(iomodel);
    559                         for(i=0;i<iomodel->numberofvertices;i++){
    560                                 if(iomodel->my_vertices[i]){
    561                                         nodes->AddObject(new Node(id0+i+1,i,lid++,0,i,false,iomodel,analysis,FSvelocityEnum));
    562                                 }
    563                         }
    564                         for(i=0;i<iomodel->numberofedges;i++){
    565                                 if(iomodel->my_edges[i]){
    566                                         nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    567                                 }
    568                         }
    569                         id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
    570                         if(iomodel->meshelementtype==PentaEnum){
    571                                 FacesPartitioning(iomodel);
    572                                 for(i=0;i<iomodel->numberoffaces;i++){
    573                                         if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
    574                                                 if(iomodel->my_faces[i]){
    575                                                         node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum);
    576                                                         nodes->AddObject(node);
    577                                                 }
    578                                         }
    579                                         else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
    580                                                 /*Nothing*/
    581                                         }
    582                                         else{
    583                                                 _error_("not supported");
    584                                         }
    585                                 }
    586                                 id0 = id0+iomodel->numberoffaces;
    587                         }
    588                         for(i=0;i<iomodel->numberofelements;i++){
    589                                 if(iomodel->my_elements[i]){
    590                                         nodes->AddObject(new Node(id0+i+1,id0+i,lid++,0,0,false,iomodel,analysis,FSvelocityEnum));
    591                                 }
    592                         }
    593 
    594                         /*No pressure*/
    595                         break;
    596 
    597                 default:
    598                         _error_("Finite element "<<EnumToStringx(finite_element)<<" not supported yet");
    599         }
    600 
    601         /*Clean up*/
    602         xDelete<bool>(my_nodes);
     674                                }
     675                        }
     676                        nodes->AddObject(node);
     677                }
     678        }
     679        xDelete<int>(nodes_approx);
     680        xDelete<int>(nodes_ranks);
     681        xDelete<int>(nodes_lids);
     682        xDelete<int>(nodes_pids);
     683
     684        nodes->common_send=common_send;
     685        nodes->common_recv=common_recv;
     686        nodes->common_send_ids=common_send_ids;
     687        nodes->common_recv_ids=common_recv_ids;
    603688}

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

    r23516 r23576  
    2222void UpdateElementsAndMaterialsDakota(Elements* elements,Materials* materials, IoModel* iomodel);
    2323void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel);
    24 void CreateNodes(Nodes*nodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum);
     24void CreateNodes(Nodes*nodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum,int* approximations=NULL);
    2525
    2626/*partitioning: */

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

    r15861 r23576  
    2020        /*Ensure that only for each cpu, the partition border nodes only will be taken into account once
    2121         * across the cluster. To do so, we flag all the clone nodes: */
    22         nodes->FlagClones(configuration_type);
     22        //nodes->FlagClones(configuration_type); /*Not needed anymore!*/
    2323
    2424        /*Go through all nodes, and build degree of freedom lists. Each node gets a fixed number of dofs. When
Note: See TracChangeset for help on using the changeset viewer.