Index: /issm/trunk-jpl/src/c/classes/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15582) +++ /issm/trunk-jpl/src/c/classes/Node.cpp (revision 15583) @@ -22,11 +22,10 @@ } /*}}}*/ -/*FUNCTION Node::Node(int node_id,int node_sid,int io_index, IoModel* iomodel,int analysis_type) {{{*/ -Node::Node(int node_id,int node_sid,int io_index, IoModel* iomodel,int analysis_type){ +/*FUNCTION Node::Node(int node_id,int node_sid,int io_index, IoModel* iomodel,int analysis_type,int approximation) {{{*/ +Node::Node(int node_id,int node_sid,int io_index, IoModel* iomodel,int analysis_type,int in_approximation){ /*Intermediary*/ int k,l; int gsize; - int node_type; /*id: */ @@ -40,10 +39,9 @@ /*indexing:*/ - node_type = reCast(iomodel->Data(FlowequationVertexEquationEnum)[io_index]); - DistributeNumDofs(&this->indexing,analysis_type,node_type); //number of dofs per node + DistributeNumDofs(&this->indexing,analysis_type,in_approximation); //number of dofs per node gsize=this->indexing.gsize; if(analysis_type==DiagnosticHorizAnalysisEnum) - this->approximation=reCast(node_type); + this->approximation=in_approximation; else this->approximation=0; @@ -55,6 +53,4 @@ if(iomodel->Data(MaskVertexongroundediceEnum)) this->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,reCast(iomodel->Data(MaskVertexongroundediceEnum)[io_index]))); - if(analysis_type==DiagnosticHorizAnalysisEnum) - this->approximation=reCast(node_type); /*set single point constraints: */ @@ -77,16 +73,16 @@ _assert_(iomodel->Data(MeshVertexonbedEnum)); _assert_(iomodel->Data(FlowequationVertexEquationEnum)); - if(node_type==SSAApproximationEnum && !reCast(iomodel->Data(MeshVertexonbedEnum)[io_index])){ + if(in_approximation==SSAApproximationEnum && !reCast(iomodel->Data(MeshVertexonbedEnum)[io_index])){ this->Deactivate(); } - if(node_type==L1L2ApproximationEnum && !reCast(iomodel->Data(MeshVertexonbedEnum)[io_index])){ + if(in_approximation==L1L2ApproximationEnum && !reCast(iomodel->Data(MeshVertexonbedEnum)[io_index])){ this->Deactivate(); } - if(node_type==SSAHOApproximationEnum && reCast(iomodel->Data(FlowequationBorderSSAEnum)[io_index])){ + if(in_approximation==SSAHOApproximationEnum && reCast(iomodel->Data(FlowequationBorderSSAEnum)[io_index])){ if(!reCast(iomodel->Data(MeshVertexonbedEnum)[io_index])){ this->Deactivate(); } } - if(node_type==SSAFSApproximationEnum && reCast(iomodel->Data(FlowequationBorderSSAEnum)[io_index])){ + if(in_approximation==SSAFSApproximationEnum && reCast(iomodel->Data(FlowequationBorderSSAEnum)[io_index])){ if(!reCast(iomodel->Data(MeshVertexonbedEnum)[io_index])){ for(k=1;k<=2;k++) this->FreezeDof(k); @@ -95,5 +91,5 @@ } /*spc all nodes on SIA*/ - if(node_type==SIAApproximationEnum){ + if(in_approximation==SIAApproximationEnum){ this->Deactivate(); } @@ -105,5 +101,5 @@ _assert_(iomodel->Data(FlowequationVertexEquationEnum)); /*Constrain all nodes that are not SIA*/ - if(reCast(node_type)!=SIAApproximationEnum){ + if(reCast(in_approximation)!=SIAApproximationEnum){ this->Deactivate(); } @@ -143,9 +139,10 @@ _printf_("Node:\n"); - _printf_(" id: " << id << "\n"); + _printf_(" id : " << id << "\n"); _printf_(" sid: " << sid << "\n"); _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); + _printf_(" approximation: " << EnumToStringx(approximation) << "\n"); indexing.Echo(); - _printf_(" inputs: " << inputs << "\n"); + _printf_(" inputs: " << inputs << "\n"); } @@ -451,4 +448,8 @@ * to a fixed value during computations. */ + if(dof>=this->indexing.gsize){ + printf("dof spc = %i\n",dof); + this->Echo(); + } _assert_(dofindexing.gsize); Index: /issm/trunk-jpl/src/c/classes/Node.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Node.h (revision 15582) +++ /issm/trunk-jpl/src/c/classes/Node.h (revision 15583) @@ -11,4 +11,5 @@ #include "../shared/shared.h" #include "./DofIndexing.h" +#include "./Update.h" class Inputs; class Hook; @@ -20,5 +21,4 @@ class ElementVector; class ElementMatrix; -#include "Update.h" /*}}}*/ @@ -38,17 +38,17 @@ IssmDouble coord_system[3][3]; - /*Node constructors, destructors {{{*/ + /*Node constructors, destructors*/ Node(); - Node(int node_id,int node_sid,int io_index, IoModel* iomodel,int analysis_type); + Node(int node_id,int node_sid,int io_index, IoModel* iomodel,int analysis_type,int approximation_in); ~Node(); - /*}}}*/ - /*Object virtual functions definitions:{{{ */ + + /*Object virtual functions definitions:*/ void Echo(); void DeepEcho(); int Id(); int ObjectEnum(); - Object *copy() {_error_("Not implemented yet (similar to Elements)"); }; - /*}}}*/ - /*Update virtual functions definitions: {{{*/ + Object *copy(){_error_("Not implemented yet (similar to Elements)"); }; + + /*Update virtual functions definitions:*/ void InputUpdateFromVector(IssmDouble* vector, int name, int type); void InputUpdateFromVector(int* vector, int name, int type); @@ -63,6 +63,6 @@ void InputUpdateFromSolution(IssmDouble* solution){_error_("Not implemented yet!");} void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("Not implemented yet!");} - /*}}}*/ - /*Node numerical routines {{{*/ + + /*Node numerical routines*/ void CreateNodalConstraints(Vector* ys); void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); @@ -97,5 +97,4 @@ void UpdateCloneDofs(int* alltruerows,int ncols,int setenum); void SetClone(int* minranks); - /*}}}*/ }; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 15582) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 15583) @@ -9,5 +9,5 @@ #include "./ModelProcessorx.h" -void CreateNodes(Nodes** pnodes, IoModel* iomodel,int analysis,int finite_element){ +void CreateNodes(Nodes** pnodes, IoModel* iomodel,int analysis,int finite_element,int approximation){ /*Intermediaries*/ @@ -26,5 +26,5 @@ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis)); + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,approximation)); } } @@ -36,5 +36,5 @@ for(j=0;j<3;j++){ if(my_nodes[3*i+j]){ - nodes->AddObject(new Node(iomodel->nodecounter+3*i+j+1,iomodel->nodecounter+3*i+j,iomodel->elements[+3*i+j]-1,iomodel,analysis)); + nodes->AddObject(new Node(iomodel->nodecounter+3*i+j+1,iomodel->nodecounter+3*i+j,iomodel->elements[+3*i+j]-1,iomodel,analysis,approximation)); } @@ -48,10 +48,10 @@ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis)); + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,approximation)); } } for(i=0;inumberofedges;i++){ if(my_edges[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis)); + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,approximation)); } } @@ -59,8 +59,9 @@ case MINIcondensedEnum: + _assert_(approximation==FSApproximationEnum); /*P1 velocity*/ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis)); + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,approximation)); } } @@ -68,5 +69,5 @@ for(i=0;inumberofvertices;i++){ if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,i,i,iomodel,analysis)); + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,i,i,iomodel,analysis,approximation)); } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15582) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15583) @@ -12,6 +12,6 @@ /*Intermediary*/ - bool isSSA,isL1L2,isHO,isFS; - int finiteelementssa; + bool isSSA,isL1L2,isHO,isFS,iscoupling; + int temp,finiteelement=-1,approximation=-1; /*Fetch parameters: */ @@ -20,22 +20,65 @@ iomodel->Constant(&isHO,FlowequationIsHOEnum); iomodel->Constant(&isFS,FlowequationIsFSEnum); - iomodel->Constant(&finiteelementssa,FlowequationFeSSAEnum); - /*Now, is the flag macayaealHO on? otherwise, do nothing: */ + /*Now, check that we have non SIA elements */ if(!isSSA & !isL1L2 & !isHO & !isFS) return; - /*Create nodes: */ - iomodel->FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, - MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); - if(finiteelementssa==0){ - CreateNodes(pnodes,iomodel,DiagnosticHorizAnalysisEnum,P1Enum); - } - else if(finiteelementssa==1){ - CreateNodes(pnodes,iomodel,DiagnosticHorizAnalysisEnum,P2Enum); + /*Do we have coupling*/ + if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.) + iscoupling = true; + else + iscoupling = false; + + /*If no coupling, call Regular CreateNodes, else, use P1 elements only*/ + if(!iscoupling){ + + /*Get finite element type*/ + if(isSSA){ + approximation=SSAApproximationEnum; + iomodel->Constant(&temp,FlowequationFeSSAEnum); + switch(temp){ + case 0 : finiteelement = P1Enum; break; + case 1 : finiteelement = P2Enum; break; + default: _error_("finite element "<FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); + CreateNodes(pnodes,iomodel,DiagnosticHorizAnalysisEnum,finiteelement,approximation); + iomodel->DeleteData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); } else{ - _error_("finite element not supported yet"); + /*Coupling: we are going to create P1 Elements only*/ + + /*First create nodes*/ + Nodes* nodes=*pnodes; + if(!nodes) nodes = new Nodes(); + + iomodel->FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,DiagnosticHorizAnalysisEnum,reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]))); + } + } + iomodel->DeleteData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); + + /*Assign output pointer: */ + *pnodes=nodes; } - iomodel->DeleteData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, - MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateNodesDiagnosticHutter.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateNodesDiagnosticHutter.cpp (revision 15582) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateNodesDiagnosticHutter.cpp (revision 15583) @@ -18,6 +18,19 @@ if(!isSIA) return; - iomodel->FetchData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); - CreateNodes(pnodes,iomodel,DiagnosticSIAAnalysisEnum,P1Enum); - iomodel->DeleteData(6,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,FlowequationVertexEquationEnum,MaskVertexonwaterEnum); + /*First create nodes*/ + Nodes* nodes=*pnodes; + if(!nodes) nodes = new Nodes(); + + iomodel->FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); + for(int i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,DiagnosticSIAAnalysisEnum,reCast(iomodel->Data(FlowequationVertexEquationEnum)[i]))); + } + } + iomodel->DeleteData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, + MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); + + /*Assign output pointer: */ + *pnodes=nodes; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 15582) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 15583) @@ -23,5 +23,5 @@ void CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel); -void CreateNodes(Nodes** pnodes, IoModel* iomodel,int analysis,int finite_element); +void CreateNodes(Nodes** pnodes, IoModel* iomodel,int analysis,int finite_element,int approximation=NoneApproximationEnum); /*Creation of fem datasets: specialised drivers: */