Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15439) @@ -947,5 +947,5 @@ /*FUNCTION Penta::GetConnectivityList {{{*/ void Penta::GetConnectivityList(int* connectivity){ - for(int i=0;iGetConnectivity(); + for(int i=0;iConnectivity(); } /*}}}*/ @@ -6308,6 +6308,6 @@ /*Find connectivity for the two nodes*/ - connectivity[0]=nodes[i]->GetConnectivity(); - connectivity[1]=nodes[i+3]->GetConnectivity(); + connectivity[0]=vertices[i]->Connectivity(); + connectivity[1]=vertices[i+3]->Connectivity(); one0=1/(IssmDouble)connectivity[0]; one1=1/(IssmDouble)connectivity[1]; @@ -7574,6 +7574,6 @@ } - connectivity[0]=nodes[node0]->GetConnectivity(); - connectivity[1]=nodes[node1]->GetConnectivity(); + connectivity[0]=vertices[node0]->Connectivity(); + connectivity[1]=vertices[node1]->Connectivity(); /*Loop on the Gauss points: */ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 15439) @@ -690,16 +690,17 @@ void Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){ - int i,j; + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = this->NumberofNodes(); + + /*First, figure out size of doflist and create it: */ + int numberofdofs=0; + for(int i=0;iGetNumberOfDofs(approximation_enum,setenum); + + /*Allocate output*/ + int* doflist=xNew(numberofdofs); + + /*Populate: */ int count=0; - int numberofdofs=0; - int* doflist=NULL; - - /*First, figure out size of doflist and create it: */ - for(i=0;i<3;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); - doflist=xNew(numberofdofs); - - /*Populate: */ - count=0; - for(i=0;i<3;i++){ + for(int i=0;iGetDofList(doflist+count,approximation_enum,setenum); count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); @@ -906,4 +907,6 @@ case P1DGEnum: return 3; + case P2Enum: + return 6; default: _error_("Element type "<element_type)<<" not supported yet"); @@ -1044,5 +1047,5 @@ /*FUNCTION Tria::GetConnectivityList {{{*/ void Tria::GetConnectivityList(int* connectivity){ - for(int i=0;iGetConnectivity(); + for(int i=0;iConnectivity(); } /*}}}*/ @@ -2182,6 +2185,5 @@ /*Intermediaries*/ - int i ,j; - int tria_node_ids[3]; + int i,j; int tria_vertex_ids[3]; int tria_type; @@ -2189,5 +2191,8 @@ IssmDouble yts; int progstabilization,balancestabilization; + int fe_ssa; bool dakota_analysis; + int numnodes; + int* tria_node_ids = NULL; /*Checks if debuging*/ @@ -2200,13 +2205,15 @@ iomodel->Constant(&progstabilization,PrognosticStabilizationEnum); iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum); + iomodel->Constant(&fe_ssa,FlowequationFeSsaEnum); iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum); /*Recover element type*/ if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){ - /*P1 Discontinuous Galerkin*/ tria_type=P1DGEnum; } + else if(analysis_type==DiagnosticHorizAnalysisEnum && fe_ssa==1){ + tria_type=P2Enum; + } else{ - /*P1 Continuous Galerkin*/ tria_type=P1Enum; } @@ -2219,19 +2226,36 @@ /*Recover nodes ids needed to initialize the node hook.*/ - if (tria_type==P1DGEnum){ - /*Discontinuous Galerkin*/ - tria_node_ids[0]=iomodel->nodecounter+3*index+1; - tria_node_ids[1]=iomodel->nodecounter+3*index+2; - tria_node_ids[2]=iomodel->nodecounter+3*index+3; - } - else{ - /*Continuous Galerkin*/ - for(i=0;i<3;i++){ - tria_node_ids[i]=iomodel->nodecounter+reCast(*(iomodel->elements+3*index+i)); //ids for vertices are in the elements array from Matlab - } + switch(tria_type){ + case P1Enum: + numnodes = 3; + tria_node_ids = xNew(numnodes); + tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0]; + tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1]; + tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2]; + break; + case P1DGEnum: + numnodes = 3; + tria_node_ids = xNew(numnodes); + tria_node_ids[0]=iomodel->nodecounter+3*index+1; + tria_node_ids[1]=iomodel->nodecounter+3*index+2; + tria_node_ids[2]=iomodel->nodecounter+3*index+3; + break; + case P2Enum: + numnodes = 6; + tria_node_ids = xNew(numnodes); + tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0]; + tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1]; + tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2]; + tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]; + tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]; + tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]; + break; + default: + _error_("Finite element "<SetHookNodes(tria_node_ids,3,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + this->SetHookNodes(tria_node_ids,numnodes,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type + xDelete(tria_node_ids); /*Fill with IoModel*/ @@ -2905,5 +2929,5 @@ /*Create Element matrix*/ for(i=0;iGetConnectivity(); + connectivity=vertices[i]->Connectivity(); Ke->values[(2*i)*numdof +(2*i) ]=1/(IssmDouble)connectivity; Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(IssmDouble)connectivity; @@ -3050,5 +3074,5 @@ gauss->GaussVertex(i); - connectivity=nodes[i]->GetConnectivity(); + connectivity=vertices[i]->Connectivity(); thickness_input->GetInputValue(&thickness,gauss); @@ -5902,5 +5926,5 @@ gauss->GaussVertex(iv); - connectivity = IssmDouble(nodes[iv]->GetConnectivity()); + connectivity = IssmDouble(vertices[iv]->Connectivity()); residual_input->GetInputValue(&residual,gauss); pe->values[iv]+=residual/connectivity; Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp (revision 15439) @@ -43,5 +43,5 @@ void TriaRef::SetElementType(int type,int type_counter){ - _assert_(type==P1Enum || type==P1DGEnum); + _assert_(type==P1Enum || type==P1DGEnum || type==P2Enum); /*initialize element type*/ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15439) @@ -288,5 +288,4 @@ /*configure elements, loads and nodes, for this new analysis: */ this->elements->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters); - this->nodes->SetCurrentConfiguration(elements,loads, nodes,vertices, materials,parameters); this->loads->SetCurrentConfiguration(elements, loads, nodes,vertices, materials,parameters); Index: /issm/trunk-jpl/src/c/classes/Hook.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Hook.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Hook.cpp (revision 15439) @@ -63,5 +63,5 @@ /*FUNCTION Hook::Echo{{{*/ void Hook::Echo(void){ - + _assert_(this); int i; if (num){ @@ -69,5 +69,5 @@ _printf_(" num=" << this->num << "\n"); _printf_(" ids: "); - for (i=0;inum;i++) _printf_(this->ids[i] << " "); + for(i=0;inum;i++) _printf_(this->ids[i] << " "); _printf_("\n"); _printf_(" offsets: "); Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15439) @@ -20,5 +20,4 @@ this->approximation=0; this->inputs=NULL; - this->hvertex=NULL; } /*}}}*/ @@ -42,7 +41,4 @@ DistributeNumDofs(&this->indexing,analysis_type,iomodel->Data(FlowequationVertexEquationEnum)+io_index); //number of dofs per node gsize=this->indexing.gsize; - - /*Hooks*/ - this->hvertex=new Hook(&vertex_id,1); //node id is the same as the vertex id, continuous galerkin! if (analysis_type==DiagnosticHorizAnalysisEnum) @@ -136,5 +132,4 @@ Node::~Node(){ delete inputs; - delete hvertex; return; } @@ -150,5 +145,4 @@ _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); indexing.Echo(); - _printf_(" hvertex: not displayed\n"); _printf_(" inputs: " << inputs << "\n"); @@ -163,6 +157,4 @@ _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); indexing.DeepEcho(); - _printf_("Vertex:\n"); - hvertex->DeepEcho(); _printf_(" inputs\n"); @@ -181,16 +173,4 @@ /*Node management:*/ -/*FUNCTION Node::Configure {{{*/ -void Node::Configure(DataSet* nodesin,Vertices* verticesin){ - - /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective - * datasets, using internal ids and off_sets hidden in hooks: */ - hvertex->configure(verticesin); - -}/*}}}*/ -/*FUNCTION Node::SetCurrentConfiguration {{{*/ -void Node::SetCurrentConfiguration(DataSet* nodesin,Vertices* verticesin){ - -}/*}}}*/ /*FUNCTION Node::GetDof {{{*/ int Node::GetDof(int dofindex,int setenum){ @@ -515,12 +495,4 @@ return approximation; -} -/*}}}*/ -/*FUNCTION Node::GetConnectivity {{{*/ -int Node::GetConnectivity(){ - - Vertex* vertex=NULL; - vertex=(Vertex*)hvertex->delivers(); - return vertex->connectivity; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Node.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.h (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Node.h (revision 15439) @@ -34,5 +34,4 @@ DofIndexing indexing; - Hook *hvertex; Inputs *inputs; //properties of this node int analysis_type; @@ -66,5 +65,4 @@ /*}}}*/ /*Node numerical routines {{{*/ - void Configure(DataSet* nodes,Vertices* vertices); void CreateNodalConstraints(Vector* ys); void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); @@ -83,5 +81,4 @@ int GetDof(int dofindex,int setenum); void CreateVecSets(Vector* pv_g,Vector* pv_f,Vector* pv_s); - int GetConnectivity(); void GetDofList(int* poutdoflist,int approximation_enum,int setenum); void GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum); Index: /issm/trunk-jpl/src/c/classes/Nodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodes.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Nodes.cpp (revision 15439) @@ -32,19 +32,4 @@ /*Numerics*/ -/*FUNCTION Nodes::Configure{{{*/ -void Nodes::Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){ - - vector::iterator object; - Node* node=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - node=dynamic_cast(*object); - node->Configure(nodes,vertices); - - } - -} -/*}}}*/ /*FUNCTION Nodes::DistributeDofs{{{*/ void Nodes::DistributeDofs(int analysis_type,int setenum){ @@ -340,16 +325,2 @@ } /*}}}*/ -/*FUNCTION Nodes::SetCurrentConfiguration{{{*/ -void Nodes::SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters){ - - vector::iterator object; - Node* node=NULL; - - for (object=objects.begin() ; object < objects.end(); object++ ){ - - node=dynamic_cast(*object); - node->SetCurrentConfiguration(nodes,vertices); - - } -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Nodes.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Nodes.h (revision 15438) +++ /issm/trunk-jpl/src/c/classes/Nodes.h (revision 15439) @@ -25,5 +25,4 @@ /*numerics*/ - void Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); void DistributeDofs(int analysis_type,int SETENUM); void FlagClones(int analysis_type); @@ -34,5 +33,4 @@ int NumberOfNodes(void); void Ranks(int* ranks,int analysis_type); - void SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); }; Index: /issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/modules/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 15439) @@ -35,12 +35,4 @@ } } - if(VerboseMProcessor()) _printf0_(" Configuring nodes...\n"); - for (i=0;iSize();i++){ - node=(Node*)nodes->GetObjectByOffset(i); - if(node->InAnalysis(configuration_type)){ - node->Configure(nodes,vertices); - } - } - if(VerboseMProcessor()) _printf0_(" Configuring materials...\n"); for (i=0;iSize();i++){ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementToEdgeConnectivity.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementToEdgeConnectivity.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementToEdgeConnectivity.cpp (revision 15439) @@ -27,17 +27,16 @@ int* element_edge_connectivity = xNew(iomodel->numberofelements*3); /*edge1 edge2 edge3*/ - /*Go through all edges and create connectivity table*/ for(int i=0;inumberofedges;i++){ /*Get the two vertices of current edge*/ - v1 = iomodel->edges[i*4+0]; _assert_(v1>=0 && v1numberofvertices); - v2 = iomodel->edges[i*4+1]; _assert_(v2>=0 && v2numberofvertices); - e1 = iomodel->edges[i*4+2]; _assert_(e1>=0 && e1numberofelements); - e2 = iomodel->edges[i*4+3]; _assert_(e2>=0 && e2numberofelements); + v1 = iomodel->edges[i*4+0]-1; _assert_(v1>=0 && v1numberofvertices); + v2 = iomodel->edges[i*4+1]-1; _assert_(v2>=0 && v2numberofvertices); + e1 = iomodel->edges[i*4+2]-1; _assert_(e1>=0 && e1numberofelements); + e2 = iomodel->edges[i*4+3]-1; _assert_(e2>=-2 && e2numberofelements); /*Process element by element*/ for(int j=0;j<3;j++){ - v3 = iomodel->elements[e1*3+j]; + v3 = iomodel->elements[e1*3+j]-1; if(v1!=v3 && v2!=v3){ element_edge_connectivity[e1*3+j]=i; @@ -45,9 +44,11 @@ } } - for(int j=0;j<3;j++){ - v3 = iomodel->elements[e2*3+j]; - if(v1!=v3 && v2!=v3){ - element_edge_connectivity[e2*3+j]=i; - break; + if(e2>-1){ + for(int j=0;j<3;j++){ + v3 = iomodel->elements[e2*3+j]-1; + if(v1!=v3 && v2!=v3){ + element_edge_connectivity[e2*3+j]=i; + break; + } } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15438) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15439) @@ -14,4 +14,5 @@ bool continuous_galerkin=true; bool isstokes,isl1l2,ismacayealpattyn; + int finiteelementssa; /*Fetch parameters: */ @@ -19,4 +20,5 @@ iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); + iomodel->Constant(&finiteelementssa,FlowequationFeSsaEnum); /*Recover pointer: */ @@ -48,8 +50,9 @@ } - if(false){ + if(finiteelementssa==1){ /*Quadratic element*/ CreateEdges(iomodel); + CreateElementToEdgeConnectivity(iomodel); int element1,element2; bool my_edge; @@ -69,5 +72,5 @@ /*Add node on edge*/ if(my_edge){ - nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,iomodel->numberofvertices+i+1,i,iomodel,DiagnosticHorizAnalysisEnum)); + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,0,iomodel,DiagnosticHorizAnalysisEnum)); } }