Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 24089) @@ -1318,4 +1318,21 @@ }/*}}}*/ +Node* Element::GetNode(int nodeindex){/*{{{*/ + _assert_(nodeindex>=0); + _assert_(nodeindexGetNumberOfNodes(this->element_type)); + return this->nodes[nodeindex]; +}/*}}}*/ +int Element::GetNodeIndex(Node* node){/*{{{*/ + + _assert_(this->nodes); + int numnodes = this->GetNumberOfNodes(this->element_type); + + for(int i=0;i* interps); + Node* GetNode(int nodeindex); + int GetNodeIndex(Node* node); void GetNodesLidList(int* lidlist); void GetNodesSidList(int* sidlist); @@ -233,9 +235,9 @@ virtual IssmDouble GetIcefrontArea(){_error_("not implemented");}; virtual void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum)=0; - virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; + virtual void GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");}; + virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");}; virtual void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level)=0; virtual void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues)=0; - virtual Node* GetNode(int node_number)=0; - virtual int GetNodeIndex(Node* node)=0; + virtual int GetVertexIndex(Vertex* vertex){_error_("not implemented");};; virtual int GetNumberOfNodes(void)=0; virtual int GetNumberOfNodes(int enum_type)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24089) @@ -1221,6 +1221,22 @@ if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria"); + int index = this->GetNodeIndex(node); + GaussPenta* gauss=new GaussPenta(); - gauss->GaussVertex(this->GetNodeIndex(node)); + gauss->GaussNode(this->element_type,index); + + input->GetInputValue(pvalue,gauss); + delete gauss; +} +/*}}}*/ +void Penta::GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){/*{{{*/ + + Input* input=inputs->GetInput(enumtype); + if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria"); + + int index = this->GetVertexIndex(vertex); + + GaussPenta* gauss=new GaussPenta(); + gauss->GaussVertex(index); input->GetInputValue(pvalue,gauss); @@ -1313,20 +1329,12 @@ } /*}}}*/ -Node* Penta::GetNode(int node_number){/*{{{*/ - _assert_(node_number>=0); - _assert_(node_numberNumberofNodes(this->element_type)); - return this->nodes[node_number]; -} -/*}}}*/ -int Penta::GetNodeIndex(Node* node){/*{{{*/ - - _assert_(nodes); - int numnodes = this->NumberofNodes(this->element_type); - - for(int i=0;i* gradient,int control_enum,int control_index,const char* data); void GetVectorFromControlInputs(Vector* gradient,int control_enum,int control_index,const char* data,int offset); + int GetVertexIndex(Vertex* vertex); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 24089) @@ -63,9 +63,7 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); void GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum); - void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");}; + void GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");}; void GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented");}; void GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues){_error_("not implemented yet");}; - Node* GetNode(int node_number){_error_("Not implemented");}; - int GetNodeIndex(Node* node){_error_("not implemented yet");}; int GetNumberOfNodes(void); int GetNumberOfNodes(int enum_type){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 24089) @@ -256,16 +256,4 @@ input->GetInputValue(pvalue,gauss); delete gauss; -} -/*}}}*/ -int Tetra::GetNodeIndex(Node* node){/*{{{*/ - - _assert_(nodes); - int numnodes = this->NumberofNodes(this->element_type); - - for(int i=0;iGetNodeIndex(node); + GaussTria* gauss=new GaussTria(); - gauss->GaussVertex(this->GetNodeIndex(node)); + gauss->GaussNode(this->element_type,index); + + input->GetInputValue(pvalue,gauss); + delete gauss; +} +/*}}}*/ +void Tria::GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){/*{{{*/ + + Input* input=inputs->GetInput(enumtype); + if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria"); + + int index = this->GetVertexIndex(vertex); + + GaussTria* gauss=new GaussTria(); + gauss->GaussVertex(index); input->GetInputValue(pvalue,gauss); @@ -2007,20 +2023,4 @@ } /*}}}*/ -Node* Tria::GetNode(int node_number){/*{{{*/ - _assert_(node_number>=0); - _assert_(node_numberNumberofNodes(this->element_type)); - return this->nodes[node_number]; - -}/*}}}*/ -int Tria::GetNodeIndex(Node* node){/*{{{*/ - - _assert_(nodes); - for(int i=0;iIsIceInElement()) return NULL; _assert_(tria->FiniteElement()==P1Enum); - int index1=tria->GetNodeIndex(nodes[0]); - int index2=tria->GetNodeIndex(nodes[1]); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); /*Intermediaries */ @@ -488,6 +488,6 @@ if(!tria->IsIceInElement()) return NULL; _assert_(tria->FiniteElement()==P1Enum); - int index1=tria->GetNodeIndex(nodes[0]); - int index2=tria->GetNodeIndex(nodes[1]); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); /*Intermediaries */ @@ -605,6 +605,6 @@ } _assert_(tria->FiniteElement()==P1Enum); - int index1=tria->GetNodeIndex(nodes[0]); - int index2=tria->GetNodeIndex(nodes[1]); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); /*Intermediaries */ Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.cpp (revision 24089) @@ -23,4 +23,5 @@ this->parameters=NULL; this->hnode=NULL; + this->hvertex=NULL; this->node=NULL; this->helement=NULL; @@ -47,4 +48,5 @@ this->hnode=new Hook(&pengrid_node_id,1); + this->hvertex=new Hook(&pengrid_node_id,1); this->helement=new Hook(&pengrid_element_id,1); @@ -52,4 +54,5 @@ this->parameters=NULL; this->node=NULL; + this->vertex=NULL; this->element=NULL; } @@ -57,4 +60,5 @@ Moulin::~Moulin(){/*{{{*/ delete hnode; + delete hvertex; delete helement; return; @@ -77,8 +81,10 @@ /*now deal with hooks and objects: */ pengrid->hnode=(Hook*)this->hnode->copy(); + pengrid->hvertex=(Hook*)this->hvertex->copy(); pengrid->helement=(Hook*)this->helement->copy(); /*corresponding fields*/ pengrid->node =(Node*)pengrid->hnode->delivers(); + pengrid->vertex=(Vertex*)pengrid->hvertex->delivers(); pengrid->element=(Element*)pengrid->helement->delivers(); @@ -91,4 +97,5 @@ _printf_(" id: " << id << "\n"); hnode->DeepEcho(); + hvertex->DeepEcho(); helement->DeepEcho(); _printf_(" parameters\n"); @@ -112,12 +119,15 @@ if(marshall_direction==MARSHALLING_BACKWARD){ this->hnode = new Hook(); + this->hvertex = new Hook(); this->helement = new Hook(); } this->hnode->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->hvertex->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); /*corresponding fields*/ node =(Node*)this->hnode->delivers(); + vertex =(Vertex*)this->hvertex->delivers(); element=(Element*)this->helement->delivers(); } @@ -135,8 +145,10 @@ * datasets, using internal ids and offsets hidden in hooks: */ hnode->configure(nodesin); + hvertex->configure(verticesin); helement->configure(elementsin); /*Get corresponding fields*/ node=(Node*)hnode->delivers(); + vertex=(Vertex*)hvertex->delivers(); element=(Element*)helement->delivers(); @@ -252,4 +264,5 @@ /*Get Element type*/ this->hnode->reset(); + this->hvertex->reset(); this->helement->reset(); Index: /issm/trunk-jpl/src/c/classes/Loads/Moulin.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Loads/Moulin.h (revision 24089) @@ -26,9 +26,11 @@ /*Hooks*/ - Hook* hnode; //hook to 1 node + Hook* hnode; //hook to 1 node + Hook* hvertex; //hook to 1 vertex Hook* helement; //hook to 1 element /*Corresponding fields*/ Node *node; + Vertex *vertex; Element *element; Index: /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Loads/Neumannflux.cpp (revision 24089) @@ -364,6 +364,6 @@ /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - int index1=tria->GetNodeIndex(nodes[0]); - int index2=tria->GetNodeIndex(nodes[1]); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); /* Start looping on the number of gaussian points: */ @@ -408,6 +408,6 @@ /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - int index1=tria->GetNodeIndex(nodes[0]); - int index2=tria->GetNodeIndex(nodes[1]); + int index1=tria->GetVertexIndex(vertices[0]); + int index2=tria->GetVertexIndex(vertices[1]); /* Start looping on the number of gaussian points: */ Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 24089) @@ -32,5 +32,5 @@ /* Intermediary */ int pos1,pos2,pos3,pos4; - int num_nodes; + int numnodes; /*numericalflux constructor data: */ @@ -50,5 +50,4 @@ if(e2==-1){ /* Boundary edge, only one element */ - num_nodes=2; numericalflux_type=BoundaryEnum; numericalflux_elem_ids[0]=e1; @@ -56,5 +55,4 @@ else{ /* internal edge: connected to 2 elements */ - num_nodes=4; numericalflux_type=InternalEnum; numericalflux_elem_ids[0]=e1; @@ -64,4 +62,6 @@ /*FIXME: hardcode element degree for now*/ this->flux_degree= P1DGEnum; + //printf("-------------- file: Numericalflux.cpp line: %i\n",__LINE__); + //this->flux_degree= P0DGEnum; /*1: Get vertices ids*/ @@ -103,4 +103,21 @@ } + switch(this->flux_degree){ + case P0DGEnum: + if(numericalflux_type==InternalEnum) numnodes = 2; + else numnodes = 1; + for(int i=0;iid = numericalflux_id; @@ -109,5 +126,5 @@ /*Hooks: */ - this->hnodes = new Hook(numericalflux_node_ids,num_nodes); + this->hnodes = new Hook(numericalflux_node_ids,numnodes); this->hvertices = new Hook(&numericalflux_vertex_ids[0],2); this->helement = new Hook(numericalflux_elem_ids,1); // take only the first element for now @@ -363,8 +380,8 @@ void Numericalflux::ResetHooks(){/*{{{*/ - this->nodes=NULL; - this->vertices=NULL; - this->element=NULL; - this->parameters=NULL; + this->nodes = NULL; + this->vertices = NULL; + this->element = NULL; + this->parameters = NULL; /*Get Element type*/ @@ -376,4 +393,5 @@ /*}}}*/ void Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ + /*Nothing to do :)*/ } Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 24089) @@ -23,6 +23,8 @@ this->parameters=NULL; this->hnodes=NULL; + this->hvertices=NULL; this->helements=NULL; this->nodes=NULL; + this->vertices=NULL; this->elements=NULL; } @@ -63,4 +65,5 @@ /*Hooks: */ this->hnodes=new Hook(riftfront_node_ids,2); + this->hvertices=new Hook(riftfront_node_ids,2); this->helements=new Hook(riftfront_elem_ids,2); @@ -88,4 +91,5 @@ this->parameters=NULL; this->nodes= NULL; + this->vertices= NULL; this->elements= NULL; @@ -95,4 +99,5 @@ this->parameters=NULL; delete hnodes; + delete hvertices; delete helements; } @@ -119,8 +124,10 @@ /*now deal with hooks and objects: */ riftfront->hnodes=(Hook*)this->hnodes->copy(); + riftfront->hvertices=(Hook*)this->hvertices->copy(); riftfront->helements=(Hook*)this->helements->copy(); /*corresponding fields*/ riftfront->nodes =(Node**)riftfront->hnodes->deliverp(); + riftfront->vertices=(Vertex**)riftfront->hvertices->deliverp(); riftfront->elements=(Element**)riftfront->helements->deliverp(); @@ -147,4 +154,5 @@ _printf_(" id: " << id << "\n"); hnodes->DeepEcho(); + hvertices->DeepEcho(); helements->DeepEcho(); _printf_(" parameters\n"); @@ -157,4 +165,5 @@ _printf_(" id: " << id << "\n"); _printf_(" hnodes: " << hnodes << "\n"); + _printf_(" hvertices: " << hvertices << "\n"); _printf_(" helements: " << helements << "\n"); _printf_(" parameters: " << parameters << "\n"); @@ -193,12 +202,15 @@ if(marshall_direction==MARSHALLING_BACKWARD){ this->hnodes = new Hook(); + this->hvertices = new Hook(); this->helements = new Hook(); } this->hnodes->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); + this->hvertices->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); this->helements->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction); /*corresponding fields*/ nodes =(Node**)this->hnodes->deliverp(); + vertices =(Vertex**)this->hvertices->deliverp(); elements =(Element**)this->helements->deliverp(); @@ -245,8 +257,10 @@ * datasets, using internal ids and offsets hidden in hooks: */ hnodes->configure(nodesin); + hvertices->configure(verticesin); helements->configure(elementsin); /*Initialize hooked fields*/ this->nodes =(Node**)hnodes->deliverp(); + this->vertices=(Vertex**)hvertices->deliverp(); this->elements=(Element**)helements->deliverp(); @@ -344,4 +358,5 @@ this->nodes=NULL; + this->vertices=NULL; this->elements=NULL; this->parameters=NULL; @@ -349,4 +364,5 @@ /*Get Element type*/ this->hnodes->reset(); + this->hvertices->reset(); this->helements->reset(); @@ -420,12 +436,8 @@ IssmDouble penalty_offset; - /*Objects: */ - Tria *tria1 = NULL; - Tria *tria2 = NULL; - /*enum of element? */ if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!"); - tria1=(Tria*)elements[0]; - tria2=(Tria*)elements[1]; + Tria* tria1=(Tria*)elements[0]; + Tria* tria2=(Tria*)elements[1]; /*Initialize Element Matrix*/ @@ -435,6 +447,6 @@ /*Get some parameters: */ this->parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum); - tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum); - tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum); + tria1->GetInputValue(&h[0],vertices[0],ThicknessEnum); + tria2->GetInputValue(&h[1],vertices[1],ThicknessEnum); if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts"); thickness=h[0]; @@ -506,12 +518,8 @@ IssmDouble pressure_water; - /*Objects: */ - Tria *tria1 = NULL; - Tria *tria2 = NULL; - /*enum of element? */ if(elements[0]->ObjectEnum()!=TriaEnum)_error_("only Tria element allowed for Riftfront load!"); - tria1=(Tria*)elements[0]; - tria2=(Tria*)elements[1]; + Tria* tria1=(Tria*)elements[0]; + Tria* tria2=(Tria*)elements[1]; /*Initialize Element Matrix*/ @@ -523,10 +531,10 @@ rho_water=tria1->FindParam(MaterialsRhoSeawaterEnum); gravity=tria1->FindParam(ConstantsGEnum); - tria1->GetInputValue(&h[0],nodes[0],ThicknessEnum); - tria2->GetInputValue(&h[1],nodes[1],ThicknessEnum); + tria1->GetInputValue(&h[0],vertices[0],ThicknessEnum); + tria2->GetInputValue(&h[1],vertices[1],ThicknessEnum); if (h[0]!=h[1])_error_("different thicknesses not supported for rift fronts"); thickness=h[0]; - tria1->GetInputValue(&b[0],nodes[0],BaseEnum); - tria2->GetInputValue(&b[1],nodes[1],BaseEnum); + tria1->GetInputValue(&b[0],vertices[0],BaseEnum); + tria2->GetInputValue(&b[1],vertices[1],BaseEnum); if (b[0]!=b[1])_error_("different beds not supported for rift fronts"); bed=b[0]; @@ -571,10 +579,10 @@ /*Ok, add contribution to first node, along the normal i==0: */ - for (j=0;j<2;j++){ + for(int j=0;j<2;j++){ pe->values[j]+=pressure*normal[j]*length; } /*Add contribution to second node, along the opposite normal: i==1 */ - for (j=0;j<2;j++){ + for(int j=0;j<2;j++){ pe->values[2+j]+= -pressure*normal[j]*length; } @@ -624,8 +632,8 @@ /*First recover velocity: */ - tria1->GetInputValue(&vx1,nodes[0],VxEnum); - tria2->GetInputValue(&vx2,nodes[1],VxEnum); - tria1->GetInputValue(&vy1,nodes[0],VyEnum); - tria2->GetInputValue(&vy2,nodes[1],VyEnum); + tria1->GetInputValue(&vx1,vertices[0],VxEnum); + tria2->GetInputValue(&vx2,vertices[1],VxEnum); + tria1->GetInputValue(&vy1,vertices[0],VyEnum); + tria2->GetInputValue(&vy2,vertices[1],VyEnum); /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 24088) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 24089) @@ -28,8 +28,10 @@ /*hooks: */ Hook* hnodes; + Hook* hvertices; Hook* helements; /*Corresponding fields*/ Node **nodes; + Vertex **vertices; Element **elements; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 24088) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 24089) @@ -497,5 +497,5 @@ /*update static arrays*/ switch(finiteelement){ - case P0Enum: + case P0Enum: case P0DGEnum: switch(iv){ case 0: coord1=1./3.; coord2=1./3.; coord3=1./3.; break;