Changeset 15439


Ignore:
Timestamp:
07/05/13 13:23:01 (12 years ago)
Author:
Mathieu Morlighem
Message:

CHG: some bug fix for quadratic elements

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

Legend:

Unmodified
Added
Removed

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

    r15438 r15439  
    947947/*FUNCTION Penta::GetConnectivityList {{{*/
    948948void  Penta::GetConnectivityList(int* connectivity){
    949         for(int i=0;i<NUMVERTICES;i++) connectivity[i]=nodes[i]->GetConnectivity();
     949        for(int i=0;i<NUMVERTICES;i++) connectivity[i]=vertices[i]->Connectivity();
    950950}
    951951/*}}}*/
    … …  
    63086308
    63096309                /*Find connectivity for the two nodes*/
    6310                 connectivity[0]=nodes[i]->GetConnectivity();
    6311                 connectivity[1]=nodes[i+3]->GetConnectivity();
     6310                connectivity[0]=vertices[i]->Connectivity();
     6311                connectivity[1]=vertices[i+3]->Connectivity();
    63126312                one0=1/(IssmDouble)connectivity[0];
    63136313                one1=1/(IssmDouble)connectivity[1];
    … …  
    75747574                }
    75757575
    7576                 connectivity[0]=nodes[node0]->GetConnectivity();
    7577                 connectivity[1]=nodes[node1]->GetConnectivity();
     7576                connectivity[0]=vertices[node0]->Connectivity();
     7577                connectivity[1]=vertices[node1]->Connectivity();
    75787578
    75797579                /*Loop on the Gauss points: */

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

    r15433 r15439  
    690690void  Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){
    691691
    692         int i,j;
     692        /*Fetch number of nodes and dof for this finite element*/
     693        int numnodes = this->NumberofNodes();
     694
     695        /*First, figure out size of doflist and create it: */
     696        int numberofdofs=0;
     697        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     698
     699        /*Allocate output*/
     700        int* doflist=xNew<int>(numberofdofs);
     701
     702        /*Populate: */
    693703        int count=0;
    694         int numberofdofs=0;
    695         int* doflist=NULL;
    696 
    697         /*First, figure out size of doflist and create it: */
    698         for(i=0;i<3;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    699         doflist=xNew<int>(numberofdofs);
    700 
    701         /*Populate: */
    702         count=0;
    703         for(i=0;i<3;i++){
     704        for(int i=0;i<numnodes;i++){
    704705                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
    705706                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    … …  
    906907                case P1DGEnum:
    907908                        return 3;
     909                case P2Enum:
     910                        return 6;
    908911                default:
    909912                        _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    … …  
    10441047/*FUNCTION Tria::GetConnectivityList {{{*/
    10451048void  Tria::GetConnectivityList(int* connectivity){
    1046         for(int i=0;i<NUMVERTICES;i++) connectivity[i]=nodes[i]->GetConnectivity();
     1049        for(int i=0;i<NUMVERTICES;i++) connectivity[i]=vertices[i]->Connectivity();
    10471050}
    10481051/*}}}*/
    … …  
    21822185
    21832186        /*Intermediaries*/
    2184         int        i                   ,j;
    2185         int        tria_node_ids[3];
     2187        int        i,j;
    21862188        int        tria_vertex_ids[3];
    21872189        int        tria_type;
    … …  
    21892191        IssmDouble yts;
    21902192        int        progstabilization,balancestabilization;
     2193        int        fe_ssa;
    21912194        bool       dakota_analysis;
     2195        int        numnodes;
     2196        int*       tria_node_ids = NULL;
    21922197
    21932198        /*Checks if debuging*/
    … …  
    22002205        iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
    22012206        iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
     2207        iomodel->Constant(&fe_ssa,FlowequationFeSsaEnum);
    22022208        iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
    22032209
    22042210        /*Recover element type*/
    22052211        if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){
    2206                 /*P1 Discontinuous Galerkin*/
    22072212                tria_type=P1DGEnum;
    22082213        }
     2214        else if(analysis_type==DiagnosticHorizAnalysisEnum && fe_ssa==1){
     2215                tria_type=P2Enum;
     2216        }
    22092217        else{
    2210                 /*P1 Continuous Galerkin*/
    22112218                tria_type=P1Enum;
    22122219        }
    … …  
    22192226
    22202227        /*Recover nodes ids needed to initialize the node hook.*/
    2221         if (tria_type==P1DGEnum){
    2222                 /*Discontinuous Galerkin*/
    2223                 tria_node_ids[0]=iomodel->nodecounter+3*index+1;
    2224                 tria_node_ids[1]=iomodel->nodecounter+3*index+2;
    2225                 tria_node_ids[2]=iomodel->nodecounter+3*index+3;
    2226         }
    2227         else{
    2228                 /*Continuous Galerkin*/
    2229                 for(i=0;i<3;i++){
    2230                         tria_node_ids[i]=iomodel->nodecounter+reCast<int,IssmDouble>(*(iomodel->elements+3*index+i)); //ids for vertices are in the elements array from Matlab
    2231                 }
     2228        switch(tria_type){
     2229                case P1Enum:
     2230                        numnodes        = 3;
     2231                        tria_node_ids   = xNew<int>(numnodes);
     2232                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     2233                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     2234                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     2235                        break;
     2236                case P1DGEnum:
     2237                        numnodes        = 3;
     2238                        tria_node_ids   = xNew<int>(numnodes);
     2239                        tria_node_ids[0]=iomodel->nodecounter+3*index+1;
     2240                        tria_node_ids[1]=iomodel->nodecounter+3*index+2;
     2241                        tria_node_ids[2]=iomodel->nodecounter+3*index+3;
     2242                        break;
     2243                case P2Enum:
     2244                        numnodes        = 6;
     2245                        tria_node_ids   = xNew<int>(numnodes);
     2246                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     2247                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     2248                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     2249                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0];
     2250                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1];
     2251                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2];
     2252                        break;
     2253                default:
     2254                        _error_("Finite element "<<EnumToStringx(tria_type)<<" not supported yet");
    22322255        }
    22332256
    22342257        /*hooks: */
    2235         this->SetHookNodes(tria_node_ids,3,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
     2258        this->SetHookNodes(tria_node_ids,numnodes,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type
     2259        xDelete<int>(tria_node_ids);
    22362260
    22372261        /*Fill with IoModel*/
    … …  
    29052929        /*Create Element matrix*/
    29062930        for(i=0;i<NUMVERTICES;i++){
    2907                 connectivity=nodes[i]->GetConnectivity();
     2931                connectivity=vertices[i]->Connectivity();
    29082932                Ke->values[(2*i)*numdof  +(2*i)  ]=1/(IssmDouble)connectivity;
    29092933                Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(IssmDouble)connectivity;
    … …  
    30503074                gauss->GaussVertex(i);
    30513075
    3052                 connectivity=nodes[i]->GetConnectivity();
     3076                connectivity=vertices[i]->Connectivity();
    30533077
    30543078                thickness_input->GetInputValue(&thickness,gauss);
    … …  
    59025926                gauss->GaussVertex(iv);
    59035927
    5904                 connectivity = IssmDouble(nodes[iv]->GetConnectivity());
     5928                connectivity = IssmDouble(vertices[iv]->Connectivity());
    59055929                residual_input->GetInputValue(&residual,gauss);
    59065930                pe->values[iv]+=residual/connectivity;

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

    r15326 r15439  
    4343void TriaRef::SetElementType(int type,int type_counter){
    4444
    45         _assert_(type==P1Enum || type==P1DGEnum);
     45        _assert_(type==P1Enum || type==P1DGEnum || type==P2Enum);
    4646
    4747        /*initialize element type*/

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

    r15429 r15439  
    288288        /*configure elements, loads and nodes, for this new analysis: */
    289289        this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
    290         this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters);
    291290        this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters);
    292291

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

    r15104 r15439  
    6363/*FUNCTION Hook::Echo{{{*/
    6464void Hook::Echo(void){
    65 
     65        _assert_(this);
    6666        int i;
    6767        if (num){
    … …  
    6969                _printf_("      num=" << this->num << "\n");
    7070                _printf_("      ids: ");
    71                 for (i=0;i<this->num;i++) _printf_(this->ids[i] << " ");
     71                for(i=0;i<this->num;i++) _printf_(this->ids[i] << " ");
    7272                _printf_("\n");
    7373                _printf_("      offsets: ");

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

    r15423 r15439  
    2020        this->approximation=0;
    2121        this->inputs=NULL;
    22         this->hvertex=NULL;
    2322}
    2423/*}}}*/
    … …  
    4241        DistributeNumDofs(&this->indexing,analysis_type,iomodel->Data(FlowequationVertexEquationEnum)+io_index); //number of dofs per node
    4342        gsize=this->indexing.gsize;
    44 
    45         /*Hooks*/
    46         this->hvertex=new Hook(&vertex_id,1); //node id is the same as the vertex id, continuous galerkin!
    4743
    4844        if (analysis_type==DiagnosticHorizAnalysisEnum)
    … …  
    136132Node::~Node(){
    137133        delete inputs;
    138         delete hvertex;
    139134        return;
    140135}
    … …  
    150145        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    151146        indexing.Echo();
    152         _printf_("   hvertex:     not displayed\n");
    153147        _printf_("   inputs:      " << inputs << "\n");
    154148
    … …  
    163157        _printf_("   analysis_type: " << EnumToStringx(analysis_type) << "\n");
    164158        indexing.DeepEcho();
    165         _printf_("Vertex:\n");
    166         hvertex->DeepEcho();
    167159        _printf_("   inputs\n");
    168160
    … …  
    181173
    182174/*Node management:*/
    183 /*FUNCTION Node::Configure {{{*/
    184 void  Node::Configure(DataSet* nodesin,Vertices* verticesin){
    185 
    186         /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
    187          * datasets, using internal ids and off_sets hidden in hooks: */
    188         hvertex->configure(verticesin);
    189 
    190 }/*}}}*/
    191 /*FUNCTION Node::SetCurrentConfiguration {{{*/
    192 void  Node::SetCurrentConfiguration(DataSet* nodesin,Vertices* verticesin){
    193 
    194 }/*}}}*/
    195175/*FUNCTION Node::GetDof {{{*/
    196176int   Node::GetDof(int dofindex,int setenum){
    … …  
    515495
    516496        return approximation;
    517 }
    518 /*}}}*/
    519 /*FUNCTION Node::GetConnectivity {{{*/
    520 int Node::GetConnectivity(){
    521 
    522         Vertex*  vertex=NULL;
    523         vertex=(Vertex*)hvertex->delivers();
    524         return vertex->connectivity;
    525497}
    526498/*}}}*/

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

    r15230 r15439  
    3434
    3535                DofIndexing  indexing;
    36                 Hook        *hvertex;
    3736                Inputs      *inputs;               //properties of this node
    3837                int          analysis_type;
    … …  
    6665                /*}}}*/
    6766                /*Node numerical routines {{{*/
    68                 void   Configure(DataSet* nodes,Vertices* vertices);
    6967                void   CreateNodalConstraints(Vector<IssmDouble>* ys);
    7068                void   SetCurrentConfiguration(DataSet* nodes,Vertices* vertices);
    … …  
    8381                int    GetDof(int dofindex,int setenum);
    8482                void   CreateVecSets(Vector<IssmDouble>* pv_g,Vector<IssmDouble>* pv_f,Vector<IssmDouble>* pv_s);
    85                 int    GetConnectivity();
    8683                void   GetDofList(int* poutdoflist,int approximation_enum,int setenum);
    8784                void   GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);

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

    r15012 r15439  
    3232
    3333/*Numerics*/
    34 /*FUNCTION Nodes::Configure{{{*/
    35 void Nodes::Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
    36 
    37         vector<Object*>::iterator object;
    38         Node* node=NULL;
    39 
    40         for ( object=objects.begin() ; object < objects.end(); object++ ){
    41 
    42                 node=dynamic_cast<Node*>(*object);
    43                 node->Configure(nodes,vertices);
    44 
    45         }
    46 
    47 }
    48 /*}}}*/
    4934/*FUNCTION Nodes::DistributeDofs{{{*/
    5035void  Nodes::DistributeDofs(int analysis_type,int setenum){
    … …  
    340325}
    341326/*}}}*/
    342 /*FUNCTION Nodes::SetCurrentConfiguration{{{*/
    343 void Nodes::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){
    344 
    345         vector<Object*>::iterator object;
    346         Node* node=NULL;
    347 
    348         for (object=objects.begin() ; object < objects.end(); object++ ){
    349 
    350                 node=dynamic_cast<Node*>(*object);
    351                 node->SetCurrentConfiguration(nodes,vertices);
    352 
    353         }
    354 }
    355 /*}}}*/

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

    r15067 r15439  
    2525
    2626                /*numerics*/
    27                 void  Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
    2827                void  DistributeDofs(int analysis_type,int SETENUM);
    2928                void  FlagClones(int analysis_type);
    … …  
    3433                int   NumberOfNodes(void);
    3534                void  Ranks(int* ranks,int analysis_type);
    36                 void  SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters);
    3735
    3836};

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

    r15104 r15439  
    3535                }
    3636        }
    37         if(VerboseMProcessor()) _printf0_("      Configuring nodes...\n");
    38         for (i=0;i<nodes->Size();i++){
    39                 node=(Node*)nodes->GetObjectByOffset(i);
    40                 if(node->InAnalysis(configuration_type)){
    41                         node->Configure(nodes,vertices);
    42                 }
    43         }
    44 
    4537        if(VerboseMProcessor()) _printf0_("      Configuring materials...\n");
    4638        for (i=0;i<materials->Size();i++){

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

    r15435 r15439  
    2727        int*  element_edge_connectivity = xNew<int>(iomodel->numberofelements*3);   /*edge1   edge2   edge3*/
    2828
    29 
    3029        /*Go through all edges and create connectivity table*/
    3130        for(int i=0;i<iomodel->numberofedges;i++){
    3231
    3332                /*Get the two vertices of current edge*/
    34                 v1 = iomodel->edges[i*4+0]; _assert_(v1>=0 && v1<iomodel->numberofvertices);
    35                 v2 = iomodel->edges[i*4+1]; _assert_(v2>=0 && v2<iomodel->numberofvertices);
    36                 e1 = iomodel->edges[i*4+2]; _assert_(e1>=0 && e1<iomodel->numberofelements);
    37                 e2 = iomodel->edges[i*4+3]; _assert_(e2>=0 && e2<iomodel->numberofelements);
     33                v1 = iomodel->edges[i*4+0]-1; _assert_(v1>=0 && v1<iomodel->numberofvertices);
     34                v2 = iomodel->edges[i*4+1]-1; _assert_(v2>=0 && v2<iomodel->numberofvertices);
     35                e1 = iomodel->edges[i*4+2]-1; _assert_(e1>=0 && e1<iomodel->numberofelements);
     36                e2 = iomodel->edges[i*4+3]-1; _assert_(e2>=-2 && e2<iomodel->numberofelements);
    3837
    3938                /*Process element by element*/
    4039                for(int j=0;j<3;j++){
    41                         v3 = iomodel->elements[e1*3+j];
     40                        v3 = iomodel->elements[e1*3+j]-1;
    4241                        if(v1!=v3 && v2!=v3){
    4342                                element_edge_connectivity[e1*3+j]=i;
    … …  
    4544                        }
    4645                }
    47                 for(int j=0;j<3;j++){
    48                         v3 = iomodel->elements[e2*3+j];
    49                         if(v1!=v3 && v2!=v3){
    50                                 element_edge_connectivity[e2*3+j]=i;
    51                                 break;
     46                if(e2>-1){
     47                        for(int j=0;j<3;j++){
     48                                v3 = iomodel->elements[e2*3+j]-1;
     49                                if(v1!=v3 && v2!=v3){
     50                                        element_edge_connectivity[e2*3+j]=i;
     51                                        break;
     52                                }
    5253                        }
    5354                }

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

    r15435 r15439  
    1414        bool   continuous_galerkin=true;
    1515        bool   isstokes,isl1l2,ismacayealpattyn;
     16        int    finiteelementssa;
    1617
    1718        /*Fetch parameters: */
    … …  
    1920        iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum);
    2021        iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum);
     22        iomodel->Constant(&finiteelementssa,FlowequationFeSsaEnum);
    2123
    2224        /*Recover pointer: */
    … …  
    4850        }
    4951
    50         if(false){
     52        if(finiteelementssa==1){
    5153
    5254                /*Quadratic element*/
    5355                CreateEdges(iomodel);
     56                CreateElementToEdgeConnectivity(iomodel);
    5457                int  element1,element2;
    5558                bool my_edge;
    … …  
    6972                        /*Add node on edge*/
    7073                        if(my_edge){
    71                                 nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,iomodel->numberofvertices+i+1,i,iomodel,DiagnosticHorizAnalysisEnum));
     74                                nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,0,iomodel,DiagnosticHorizAnalysisEnum));
    7275                        }
    7376                }
Note: See TracChangeset for help on using the changeset viewer.