Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 847) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 848) @@ -1299,9 +1299,10 @@ } - -void DataSet::VelocityExtrude(Vec ug,double* ug_serial){ + +void DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){ vector::iterator object; Penta* penta=NULL; + Node* node=NULL; for ( object=objects.begin() ; object < objects.end(); object++ ){ @@ -1310,43 +1311,18 @@ penta=(Penta*)(*object); - penta->VelocityExtrude(ug,ug_serial); - - } - } - -} -void DataSet::ThicknessExtrude(Vec tg,double* tg_serial){ - - vector::iterator object; - Penta* penta=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PentaEnum()){ - - penta=(Penta*)(*object); - penta->ThicknessExtrude(tg,tg_serial); - - } - } - -} -void DataSet::VelocityExtrudeAllElements(Vec ug,double* ug_serial){ - - vector::iterator object; - Penta* penta=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PentaEnum()){ - - penta=(Penta*)(*object); - penta->VelocityExtrudeAllElements(ug,ug_serial); - - } - } - -} -void DataSet::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){ + penta->FieldExtrude(field,field_serial,field_name,collapse); + + } + if((*object)->Enum()==NodeEnum()){ + + node=(Node*)(*object); + node->FieldExtrude(field,field_serial,field_name); + + } + } + +} + +void DataSet::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){ vector::iterator object; @@ -1358,21 +1334,5 @@ node=(Node*)(*object); - node->VelocityDepthAverageAtBase(ug,ug_serial); - } - } - -} -void DataSet::SlopeExtrude(Vec sg,double* sg_serial){ - - vector::iterator object; - Penta* penta=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - if((*object)->Enum()==PentaEnum()){ - - penta=(Penta*)(*object); - penta->SlopeExtrude(sg,sg_serial); - + node->FieldDepthAverageAtBase(field,field_serial,fieldname); } } @@ -1397,2 +1357,19 @@ +void DataSet::UpdateNodePositions(double* thickness,double* bed){ + + vector::iterator object; + Node* node=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==NodeEnum()){ + + node=(Node*)(*object); + node->UpdateNodePosition(thickness,bed); + } + } + + + +} Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 847) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 848) @@ -76,12 +76,10 @@ void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type); void Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type); - void VelocityExtrude(Vec ug,double* ug_serial); - void ThicknessExtrude(Vec ug,double* ug_serial); - void VelocityExtrudeAllElements(Vec ug,double* ug_serial); - void VelocityDepthAverageAtBase(Vec ug,double* ug_serial); - void SlopeExtrude(Vec sg,double* sg_serial); + void FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname); int DeleteObject(Object* object); void ComputePressure(Vec p_g); int FindResult(void* pvalue, char* name); + void FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse); + void UpdateNodePositions(double* thickness,double* bed); }; Index: /issm/trunk/src/c/FieldDepthAveragex/FieldDepthAveragex.cpp =================================================================== --- /issm/trunk/src/c/FieldDepthAveragex/FieldDepthAveragex.cpp (revision 848) +++ /issm/trunk/src/c/FieldDepthAveragex/FieldDepthAveragex.cpp (revision 848) @@ -0,0 +1,50 @@ +/*!\file FieldDepthAveragex + * \brief: average field throfieldh thickness + */ + +#include "./FieldDepthAveragex.h" + +#undef __FUNCT__ +#define __FUNCT__ "FieldDepthAveragex" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void FieldDepthAveragex( Vec field, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,char* fieldname){ + + double* field_serial=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes, materials); + nodes->Configure(elements,loads, nodes, materials); + + /*Serialize field: */ + VecToMPISerial(&field_serial,field); + + /*Depth average field, onto base of mesh: */ + nodes->FieldDepthAverageAtBase(field,field_serial,fieldname); + + /*Assemble vector: */ + VecAssemblyBegin(field); + VecAssemblyEnd(field); + + /*ok, field is now correct at the base. We need to extrude the value of the base + * through thickness :*/ + + /*Serialize field: */ + xfree((void**)&field_serial); + VecToMPISerial(&field_serial,field); + + /*Extrude field to all nodes: */ + nodes->FieldExtrude(field,field_serial,fieldname,0); + + /*Assemble vector: */ + VecAssemblyBegin(field); + VecAssemblyEnd(field); + + /*Free ressources:*/ + xfree((void**)&field_serial); + +} Index: /issm/trunk/src/c/FieldDepthAveragex/FieldDepthAveragex.h =================================================================== --- /issm/trunk/src/c/FieldDepthAveragex/FieldDepthAveragex.h (revision 848) +++ /issm/trunk/src/c/FieldDepthAveragex/FieldDepthAveragex.h (revision 848) @@ -0,0 +1,14 @@ +/*!\file: FieldDepthAveragex.h + * \brief header file for averaging a field throfieldh thickness + */ + +#ifndef _FIELDDEPTHAVERAGEX_H +#define _FIELDDEPTHAVERAGEX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void FieldDepthAveragex( Vec field, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,char* fieldname); + +#endif /* _FIELDDEPTHAVERAGEX_H */ + Index: /issm/trunk/src/c/FieldExtrudex/FieldExtrudex.cpp =================================================================== --- /issm/trunk/src/c/FieldExtrudex/FieldExtrudex.cpp (revision 848) +++ /issm/trunk/src/c/FieldExtrudex/FieldExtrudex.cpp (revision 848) @@ -0,0 +1,46 @@ +/*!\file FieldExtrudex + * \brief: vertical velocity extrusion + */ + +#include "./FieldExtrudex.h" + +#undef __FUNCT__ +#define __FUNCT__ "FieldExtrudex" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void FieldExtrudex( Vec field, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,char* field_name,int collapse){ + + double* field_serial=NULL; + + /*First, get elements and nodes configured: */ + elements->Configure(elements,loads, nodes, materials); + nodes->Configure(elements,loads, nodes, materials); + + /*Serialize field: */ + VecToMPISerial(&field_serial,field); + + /*ok, is the extrusion dependent upon elements or just nodes? If collapse is 1, it means we only want + * to extrude the field on elemenets that are collapsed (element operation). If collapse is 0, it means + * we want to extrude field on every node (node operation): */ + + if(collapse){ + /*Extrude field vertically: */ + elements->FieldExtrude(field,field_serial,field_name,collapse); + } + else{ + /*Extrude field vertically: */ + nodes->FieldExtrude(field,field_serial,field_name,collapse); + } + + /*Assemble vector: */ + VecAssemblyBegin(field); + VecAssemblyEnd(field); + + /*Free ressources:*/ + xfree((void**)&field_serial); + +} Index: /issm/trunk/src/c/FieldExtrudex/FieldExtrudex.h =================================================================== --- /issm/trunk/src/c/FieldExtrudex/FieldExtrudex.h (revision 848) +++ /issm/trunk/src/c/FieldExtrudex/FieldExtrudex.h (revision 848) @@ -0,0 +1,14 @@ +/*!\file: FieldExtrudex.h + * \brief header file for field extrusion + */ + +#ifndef _FIELDEXTRUDEX_H +#define _FIELDEXTRUDEX_H + +#include "../DataSet/DataSet.h" + +/* local prototypes: */ +void FieldExtrudex( Vec field, DataSet* elements,DataSet* nodes, DataSet* loads, DataSet* materials,char* field_name, int collapse); + +#endif /* _FIELDEXTRUDEX_H */ + Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 847) +++ /issm/trunk/src/c/Makefile.am (revision 848) @@ -264,12 +264,8 @@ ./ProcessParamsx/ProcessParamsx.cpp\ ./ProcessParamsx/ProcessParamsx.h\ - ./ThicknessExtrudex/ThicknessExtrudex.cpp\ - ./ThicknessExtrudex/ThicknessExtrudex.h\ - ./VelocityExtrudex/VelocityExtrudex.cpp\ - ./VelocityExtrudex/VelocityExtrudex.h\ - ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\ - ./VelocityDepthAveragex/VelocityDepthAveragex.h\ - ./SlopeExtrudex/SlopeExtrudex.cpp\ - ./SlopeExtrudex/SlopeExtrudex.h + ./FieldDepthAveragex/FieldDepthAveragex.cpp\ + ./FieldDepthAveragex/FieldDepthAveragex.h\ + ./FieldExtrudex/FieldExtrudex.cpp\ + ./FieldExtrudex/FieldExtrudex.h @@ -532,16 +528,12 @@ ./Solverx/Solverx.cpp\ ./Solverx/Solverx.h\ - ./ThicknessExtrudex/ThicknessExtrudex.cpp\ - ./ThicknessExtrudex/ThicknessExtrudex.h\ ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\ ./ProcessParamsx/ProcessParamsx.cpp\ ./ProcessParamsx/ProcessParamsx.h\ - ./VelocityExtrudex/VelocityExtrudex.cpp\ - ./VelocityExtrudex/VelocityExtrudex.h\ - ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\ - ./VelocityDepthAveragex/VelocityDepthAveragex.h\ - ./SlopeExtrudex/SlopeExtrudex.cpp\ - ./SlopeExtrudex/SlopeExtrudex.h\ + ./FieldDepthAveragex/FieldDepthAveragex.cpp\ + ./FieldDepthAveragex/FieldDepthAveragex.h\ + ./FieldExtrudex/FieldExtrudex.cpp\ + ./FieldExtrudex/FieldExtrudex.h\ ./parallel/diagnostic_core.cpp\ ./parallel/diagnostic_core_linear.cpp\ Index: /issm/trunk/src/c/UpdateNodePositionsx/UpdateNodePositionsx.cpp =================================================================== --- /issm/trunk/src/c/UpdateNodePositionsx/UpdateNodePositionsx.cpp (revision 848) +++ /issm/trunk/src/c/UpdateNodePositionsx/UpdateNodePositionsx.cpp (revision 848) @@ -0,0 +1,36 @@ +/*!\file UpdateNodePositionsx + * \brief update node positions using new geometry defined by new thickness and new bed. + */ + +#include "./UpdateNodePositionsx.h" + +#undef __FUNCT__ +#define __FUNCT__ "UpdateNodePositionsx" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" + +void UpdateNodePositionsx( DataSet* elements,DataSet* nodes,DataSet* loads, DataSet* materials, Vec thickness,Vec bed){ + + /*intermediary: */ + double* thickness_serial=NULL; + double* bed_serial=NULL; + + /*First, get elements and loads configured: */ + elements->Configure(elements,loads, nodes, materials); + loads->Configure(elements, loads, nodes, materials); + nodes->Configure(elements, loads, nodes, materials); + + /*serialize inputs: */ + VecToMPISerial(&thickness_serial,thickness); + VecToMPISerial(&bed_serial,bed); + + /*Update elements, nodes, loads and materials from new geometry: */ + nodes->UpdateNodePositions(thickness_serial,bed_serial); + + /*Free ressources:*/ + xfree((void**)&thickness_serial); + xfree((void**)&bed_serial); +} Index: /issm/trunk/src/c/UpdateNodePositionsx/UpdateNodePositionsx.h =================================================================== --- /issm/trunk/src/c/UpdateNodePositionsx/UpdateNodePositionsx.h (revision 848) +++ /issm/trunk/src/c/UpdateNodePositionsx/UpdateNodePositionsx.h (revision 848) @@ -0,0 +1,15 @@ +/*!\file: UpdateNodePositionsx.h + * \brief header file for updating node positions using new geometry defined by new thickness and new bed. + */ + +#ifndef _UPDATENODEPOSITIONSXX_H +#define _UPDATENODEPOSITIONSXX_H + +#include "../DataSet/DataSet.h" +#include "../objects/ParameterInputs.h" + +/* local prototypes: */ +void UpdateNodePositionsx( DataSet* elements,DataSet* nodes,DataSet* loads, DataSet* materials, Vec thickness,Vec bed); + +#endif /* _UPDATENODEPOSITIONSXX_H */ + Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 847) +++ /issm/trunk/src/c/issm.h (revision 848) @@ -34,4 +34,5 @@ #include "./UpdateFromInputsx/UpdateFromInputsx.h" #include "./UpdateGeometryx/UpdateGeometryx.h" +#include "./UpdateNodePositionsx/UpdateNodePositionsx.h" #include "./PenaltySystemMatricesx/PenaltySystemMatricesx.h" #include "./Reducematrixfromgtofx/Reducematrixfromgtofx.h" @@ -45,10 +46,8 @@ #include "./Misfitx/Misfitx.h" #include "./ControlConstrainx/ControlConstrainx.h" -#include "./VelocityExtrudex/VelocityExtrudex.h" -#include "./VelocityDepthAveragex/VelocityDepthAveragex.h" +#include "./FieldDepthAveragex/FieldDepthAveragex.h" #include "./GriddataMeshToGridx/GriddataMeshToGridx.h" #include "./ComputePressurex/ComputePressurex.h" -#include "./SlopeExtrudex/SlopeExtrudex.h" -#include "./ThicknessExtrudex/ThicknessExtrudex.h" +#include "./FieldExtrudex/FieldExtrudex.h" Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 847) +++ /issm/trunk/src/c/objects/Node.cpp (revision 848) @@ -550,6 +550,6 @@ #undef __FUNCT__ -#define __FUNCT__ "Node::VelocityDepthAverageAtBase" -void Node::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){ +#define __FUNCT__ "Node::FieldDepthAverageAtBase" +void Node::FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname){ /* node data: */ @@ -560,11 +560,8 @@ Node* node=NULL; Node* upper_node=NULL; - double velocity2[2]; - double velocity1[2]; - double velocity_average[2]; - double sum[2]; double z1,z2,dz; double thickness; + /*Are we on the base, not on the surface, and not on a clone node?:*/ if(onbed==1 & clone==0 &onsurface==0){ @@ -572,53 +569,206 @@ doflist1=GetDofList1(); - /*this node is on the bed. We are going to: - * follow the upper nodes until we reach the surface. At each upper node, we'll grab the - * velocity for this node, and add it to ug. We'll finish by - * we grab the velocity. Once we know the velocity, we follow the upper nodes, - * inserting the same velocity value into ug, until we reach the surface: */ - sum[0]=0; - sum[1]=0; - thickness=0; - - /*get dofs for this base node velocity: we know there are two dofs in ug_serial */ - dofx=2*doflist1; - dofy=2*doflist1+1; - - node=this; - for(;;){ - - if (node->IsOnSurface())break; - - doflist1=node->GetDofList1(); + /*this node is on the bed. We are going to, follow the upper nodes until we reach the surface. At each upper node, + * we'll grab the * field for this node, and add it to overall field: */ + + if(strcmp(fieldname,"velocity")==0){ + + /*field is a velocity, 2 dofs per node: */ + double velocity2[2]; + double velocity1[2]; + double velocity_average[2]; + double sum[2]; + + sum[0]=0; + sum[1]=0; + thickness=0; + + /*get dofs for this base node velocity: we know there are two dofs in field_serial */ + dofx=2*doflist1; + dofy=2*doflist1+1; + + node=this; + for(;;){ + + if (node->IsOnSurface())break; + + doflist1=node->GetDofList1(); + + velocity1[0]=field_serial[2*doflist1]; + velocity1[1]=field_serial[2*doflist1+1]; + z1=node->GetZ(); + + upper_node=node->GetUpperNode(); + doflist1=upper_node->GetDofList1(); - velocity1[0]=ug_serial[2*doflist1]; - velocity1[1]=ug_serial[2*doflist1+1]; - z1=node->GetZ(); - - upper_node=node->GetUpperNode(); - doflist1=upper_node->GetDofList1(); - - velocity2[0]=ug_serial[2*doflist1]; - velocity2[1]=ug_serial[2*doflist1+1]; - z2=upper_node->GetZ(); - - dz=(z2-z1); - thickness+=dz; - velocity_average[0]=(velocity1[0]+velocity2[0])/2.0; - velocity_average[1]=(velocity1[1]+velocity2[1])/2.0; - - sum[0]+=velocity_average[0]*dz; - sum[1]+=velocity_average[1]*dz; - - /* get next node: */ - node=node->GetUpperNode(); + velocity2[0]=field_serial[2*doflist1]; + velocity2[1]=field_serial[2*doflist1+1]; + z2=upper_node->GetZ(); + + dz=(z2-z1); + thickness+=dz; + velocity_average[0]=(velocity1[0]+velocity2[0])/2.0; + velocity_average[1]=(velocity1[1]+velocity2[1])/2.0; + + sum[0]+=velocity_average[0]*dz; + sum[1]+=velocity_average[1]*dz; + + /* get next node: */ + node=node->GetUpperNode(); + } + + sum[0]=sum[0]/thickness; + sum[1]=sum[1]/thickness; + + /* Plfield velocity_average*deltaH/H into base of field: */ + VecSetValues(field,1,&dofx,&sum[0],INSERT_VALUES); + VecSetValues(field,1,&dofy,&sum[1],INSERT_VALUES); } - - sum[0]=sum[0]/thickness; - sum[1]=sum[1]/thickness; - - /* Plug velocity_average*deltaH/H into base of ug: */ - VecSetValues(ug,1,&dofx,&sum[0],INSERT_VALUES); - VecSetValues(ug,1,&dofy,&sum[1],INSERT_VALUES); - } -} + else{ + /*field is regular, 1 dof per node: */ + double field2; + double field1; + double field_average; + double sum; + + sum=0; + thickness=0; + + /*get dofs for this base node velocity: we know there are two dofs in field_serial */ + dofx=doflist1; + + node=this; + for(;;){ + + if (node->IsOnSurface())break; + + doflist1=node->GetDofList1(); + + field1=field_serial[doflist1]; + z1=node->GetZ(); + + upper_node=node->GetUpperNode(); + doflist1=upper_node->GetDofList1(); + + field2=field_serial[2*doflist1]; + z2=upper_node->GetZ(); + + dz=(z2-z1); + thickness+=dz; + field_average=(field1+field2)/2.0; + + sum+=field_average*dz; + + /* get next node: */ + node=node->GetUpperNode(); + } + + sum=sum/thickness; + + /* Plug field_average*deltH/H into base of field: */ + VecSetValues(field,1,&dofx,&sum,INSERT_VALUES); + } + } +} + +#undef __FUNCT__ +#define __FUNCT__ "Node::FieldExtrude" +void Node::FieldExtrude(Vec field,double* field_serial,char* field_name){ + + /* node data: */ + int numberofdofspernode; + Node* node=NULL; + int i; + + /*Is this node on bed? :*/ + if (onbed){ + + if (strcmp(field_name,"velocity")==0){ + + /* node data: */ + const int numdof=2; + int doflist[numdof]; + int nodedofs[2]; + double fieldnode[2]; + + + GetDofList(&doflist[0],&numberofdofspernode); + + //initilaize node + node=this; + + /*get field for this base node: */ + fieldnode[0]=field_serial[doflist[0]]; + fieldnode[1]=field_serial[doflist[1]]; + + //go through all nodes which sit on top of this node, until we reach the surface, + //and plfield field in field + for(;;){ + + node->GetDofList(&nodedofs[0],&numberofdofspernode); + VecSetValues(field,1,&nodedofs[0],&fieldnode[0],INSERT_VALUES); + VecSetValues(field,1,&nodedofs[1],&fieldnode[1],INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } //if (strcmp(field_name,"velocity")==0) + else if ( + (strcmp(field_name,"thickness")==0) || + (strcmp(field_name,"surface")==0) || + (strcmp(field_name,"bed")==0) || + (strcmp(field_name,"slopex")==0) || + (strcmp(field_name,"slopey")==0) + ){ + + /* node data: */ + int doflist; + int nodedofs; + double fieldnode; + + GetDofList(&doflist,&numberofdofspernode); + + //initilaize node + node=this; + + /*get field for this bed node: */ + fieldnode=field_serial[doflist]; + + //go through all nodes which sit on top of this node, until we reach the surface, + //and pltg fieldnode in field: + for(;;){ + + node->GetDofList(&nodedofs,&numberofdofspernode); + VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!")); + + } //if (extrude) +} + +void Node::UpdateNodePosition(double* thickness,double* bed){ + + int dof; + + dof=this->GetDofList1(); + + if(onbed){ + + /*this node is on the bed, set its z position to the new bed: */ + this->x[2]=bed[dof]; + } + else if (onsurface){ + + /*this node is on the surface, set its z position to the new bed+ new thickness: */ + this->x[2]=bed[dof]+thickness[dof]; + } + else{ + /*do nothing, we only update grids belonging to the exterior layers: */ + } +} Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 847) +++ /issm/trunk/src/c/objects/Node.h (revision 848) @@ -81,7 +81,9 @@ int IsOnSurface(); void FreezeDof(int dof); - void VelocityDepthAverageAtBase(Vec ug,double* ug_serial); + void FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname); int IsOnShelf(); int IsOnSheet(); + void FieldExtrude(Vec field,double* field_serial,char* field_name); + void UpdateNodePosition(double* thickness,double* bed); }; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 847) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 848) @@ -1922,200 +1922,112 @@ #undef __FUNCT__ -#define __FUNCT__ "Penta::VelocityExtrude" -void Penta::VelocityExtrude(Vec ug,double* ug_serial){ - +#define __FUNCT__ "Penta::FieldExtrude" +void Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){ + /* node data: */ const int numgrids=6; - const int numdof=2*numgrids; - int doflist[numdof]; - int nodedofs[2]; int numberofdofspernode; - Node* node=NULL; int i; - double velocity[2]; - - - if((collapse==1) && (onbed==1)){ - - GetDofList(&doflist[0],&numberofdofspernode); - - /*this penta is a collapsed macayeal. For each node on the base of this penta, - * we grab the velocity. Once we know the velocity, we follow the upper nodes, - * inserting the same velocity value into ug, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get velocity for this base node: */ - velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]]; - velocity[1]=ug_serial[doflist[numberofdofspernode*i+1]]; - - //go througn all nodes which sit on top of this node, until we reach the surface, - //and plug velocity in ug - for(;;){ - - node->GetDofList(&nodedofs[0],&numberofdofspernode); - VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES); - VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); + int extrude=0; + + /*Figure out if we should extrude for this element: */ + if (iscollapsed){ + /*From higher level, we are told to extrude only elements that have the collapse flag on: */ + if (collapse)extrude=1; + else extrude=0; + } + else{ + /*From higher level, we are told to extrude all elements: */ + extrude=1; + } + + /*Now, extrusion starts from the bed on, so double check this element is on + * the bedrock: */ + if(onbed==0)extrude=0; + + /*Go on and extrude field: */ + if (extrude){ + + if (strcmp(field_name,"velocity")==0){ + + /* node data: */ + const int numdof=2*numgrids; + int doflist[numdof]; + int nodedofs[2]; + double fieldel[2]; + + + GetDofList(&doflist[0],&numberofdofspernode); + + /*this penta is a collapsed macayeal. For each node on the base of this penta, + * we grab the field. Once we know the field, we follow the upper nodes, + * inserting the same field value into field, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + + /*get field for this base node: */ + fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]]; + fieldel[1]=field_serial[doflist[numberofdofspernode*i+1]]; + + //go throfieldn all nodes which sit on top of this node, until we reach the surface, + //and plfield field in field + for(;;){ + + node->GetDofList(&nodedofs[0],&numberofdofspernode); + VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES); + VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } } + } //if (strcmp(field_name,"velocity")==0) + else if ( + (strcmp(field_name,"thickness")==0) || + (strcmp(field_name,"surface")==0) || + (strcmp(field_name,"bed")==0) || + (strcmp(field_name,"slopex")==0) || + (strcmp(field_name,"slopey")==0) + ){ + + /* node data: */ + const int numdof=1*numgrids; + int doflist[numdof]; + int nodedofs; + double fieldel; + + GetDofList(&doflist[0],&numberofdofspernode); + + /*this penta is on the bed. For each node on the base of this penta, + * we grab the thickness. Once we know the thickness, we follow the upper nodes, + * inserting the same thickness value into tg, until we reach the surface: */ + for(i=0;i<3;i++){ + + node=nodes[i]; //base nodes + + /*get velocity for this base node: */ + fieldel=field_serial[doflist[numberofdofspernode*i+0]]; + + //go through all nodes which sit on top of this node, until we reach the surface, + //and pltg fieldel in field: + for(;;){ + + node->GetDofList(&nodedofs,&numberofdofspernode); + VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES); + + if (node->IsOnSurface())break; + /*get next node: */ + node=node->GetUpperNode(); + } + } + } - - } - -} - -#undef __FUNCT__ -#define __FUNCT__ "Penta::VelocityExtrudeAllElements" -void Penta::VelocityExtrudeAllElements(Vec ug,double* ug_serial){ - - /* node data: */ - const int numgrids=6; - int doflist[numgrids]; - int doflist1; - int nodedofs[2]; - - Node* node=NULL; - int i; - double velocity[2]; - - - if(onbed==1){ - - GetDofList1(&doflist[0]); - - /*this penta is on the bed. For each node on the base of this penta, - * we grab the velocity. Once we know the velocity, we follow the upper nodes, - * inserting the same velocity value into ug, until we reach the surface: - * the velocity has two ddls vx and vy*/ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get velocity for this base node: */ - velocity[0]=ug_serial[2*doflist[i]+0]; - velocity[1]=ug_serial[2*doflist[i]+1]; - - //go througn all nodes which sit on top of this node, until we reach the surface, - //and plug velocity in ug - for(;;){ - - doflist1=node->GetDofList1(); - nodedofs[0]=2*doflist1; - nodedofs[1]=2*doflist1+1; - VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES); - VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - - } - -} - -#undef __FUNCT__ -#define __FUNCT__ "Penta::ThicknessExtrude" -void Penta::ThicknessExtrude(Vec tg,double* tg_serial){ - - /* node data: */ - const int numgrids=6; - const int numdof=1*numgrids; - int doflist[numdof]; - int nodedofs; - int numberofdofspernode; - - Node* node=NULL; - int i; - double thickness; - - - if(onbed==1){ - - GetDofList(&doflist[0],&numberofdofspernode); - - /*this penta is on the bed. For each node on the base of this penta, - * we grab the thickness. Once we know the thickness, we follow the upper nodes, - * inserting the same thickness value into tg, until we reach the surface: */ - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get velocity for this base node: */ - thickness=tg_serial[doflist[numberofdofspernode*i+0]]; - - //go through all nodes which sit on top of this node, until we reach the surface, - //and pltg velocity in tg - for(;;){ - - node->GetDofList(&nodedofs,&numberofdofspernode); - VecSetValues(tg,1,&nodedofs,&thickness,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - - } - -} - -#undef __FUNCT__ -#define __FUNCT__ "Penta::SlopeExtrude" -void Penta::SlopeExtrude(Vec sg,double* sg_serial){ - - /* node data: */ - const int numgrids=6; - const int NDOF1=1; - const int numdof=NDOF1*numgrids; - int doflist[numdof]; - int nodedof; - int numberofdofspernode; - - Node* node=NULL; - int i; - double slope; - - - if(onbed==1){ - - GetDofList(&doflist[0],&numberofdofspernode); - - if(numberofdofspernode!=1)throw ErrorException(__FUNCT__," slope can only be extruded on 1 dof per node"); - - /*For each node on the base of this penta, we grab the slope. Once we know the slope, we follow the upper nodes, - * inserting the same slope value into sg, until we reach the surface: */ - - for(i=0;i<3;i++){ - - node=nodes[i]; //base nodes - - /*get velocity for this base node: */ - slope=sg_serial[doflist[i]]; - - //go throsgn all nodes which sit on top of this node, until we reach the surface, - //and plsg slope in sg - for(;;){ - - node->GetDofList(&nodedof,&numberofdofspernode); - VecSetValues(sg,1,&nodedof,&slope,INSERT_VALUES); - - if (node->IsOnSurface())break; - /*get next node: */ - node=node->GetUpperNode(); - } - } - - } - -} - + else throw ErrorException(__FUNCT__,exprintf("%s%s%s"," field ",field_name," not supported yet!")); + + } //if (extrude) +} #undef __FUNCT__ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 847) +++ /issm/trunk/src/c/objects/Penta.h (revision 848) @@ -111,8 +111,5 @@ void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4); void GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4); - void VelocityExtrudeAllElements(Vec ug,double* ug_serial); - void VelocityExtrude(Vec ug,double* ug_serial); - void ThicknessExtrude(Vec ug,double* ug_serial); - void SlopeExtrude(Vec sg,double* sg_serial); + void FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed); void ComputePressure(Vec p_g); void CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/parallel/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 847) +++ /issm/trunk/src/c/parallel/diagnostic_core.cpp (revision 848) @@ -85,6 +85,6 @@ if(debug)_printf_("%s\n","extruding slopes in 3d..."); - SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); - SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); + FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0); + FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0); } @@ -122,5 +122,5 @@ if(debug)_printf_("%s\n"," extruding horizontal velocities..."); - VecDuplicatePatch(&ug_horiz,ug); VelocityExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials); + VecDuplicatePatch(&ug_horiz,ug); FieldExtrudex( ug_horiz,fem_dh->elements,fem_dh->nodes, fem_dh->loads,fem_dh-> materials,"velocity",1); if(debug)_printf_("%s\n"," computing vertical velocities..."); @@ -145,6 +145,6 @@ diagnostic_core_linear(&slopex,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum()); diagnostic_core_linear(&slopey,fem_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum()); - SlopeExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); - SlopeExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials); + FieldExtrudex( slopex, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopex",0); + FieldExtrudex( slopey, fem_sl->elements,fem_sl->nodes,fem_sl->loads,fem_sl->materials,"slopey",0); inputs->Add("bedslopex",slopex,numberofdofspernode_sl,numberofnodes); Index: /issm/trunk/src/c/parallel/prognostic_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic_core.cpp (revision 847) +++ /issm/trunk/src/c/parallel/prognostic_core.cpp (revision 848) @@ -39,5 +39,5 @@ _printf_("depth averaging velocity..."); u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity - VelocityDepthAveragex( u_g, fem->elements,fem->nodes, fem->loads, fem->materials); + FieldDepthAveragex( u_g, fem->elements,fem->nodes, fem->loads, fem->materials,"velocity"); inputs->Add("velocity_average",u_g,2,numberofnodes); @@ -46,5 +46,5 @@ _printf_("extrude computed thickness on all layers:\n"); - ThicknessExtrudex( h_g, fem->elements,fem->nodes, fem->loads, fem->materials); + FieldExtrudex( h_g, fem->elements,fem->nodes, fem->loads, fem->materials,"thickness",0); /*Plug results into output dataset: */