Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3868) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3869) @@ -945,4 +945,104 @@ } +} +/*}}}*/ +/*FUNCTION DataSet::DepthAverageInputAtBase(int enum_type) {{{1*/ +void DataSet::DepthAverageInputAtBase(int enum_type){ + + vector::iterator object; + Element* element=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + element=(Element*)(*object); + element->DepthAverageInputAtBase(enum_type); + } + } + +} +/*}}}*/ +/*FUNCTION DataSet::InputExtrude(int enum_type) {{{1*/ +void DataSet::InputExtrude(int enum_type){ + + vector::iterator object; + Penta* penta=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if((*object)->Enum()==PentaEnum){ + penta=(Penta*)(*object); + penta->InputExtrude(enum_type); + } + } +} +/*}}}*/ +/*FUNCTION DataSet::InputToResult{{{1*/ +void DataSet::InputToResult(Mat* psolution,int enum_type){ + + /*Output*/ + Mat solution=NULL; + + /*Intermediary*/ + int i; + vector::iterator object; + extern int num_procs; + extern int my_rank; + Element* element=NULL; + int elemcols; + int cols=0,lines=0; + int local_cols=0; + int* all_cols=NULL; + int local_lines=0; + int* all_lines=NULL; + + ISSMERROR("not implemented yet"); + /*First: Get number of columns to be allocated*/ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + if(EnumIsElement((*object)->Enum())){ + + /*increase number of lines by 1*/ + local_lines++; + + /*Get number of columns required for this element*/ + element=(Element*)(*object); + //elemcols=element->InputToResultAllocation(enum_type); + + /*Increase number of columns if higher than previous value*/ + //if(elemcols>local_cols) local_cols=elemcols; + } + } + + /*OK, sum contributions of all cpus*/ + all_lines=(int*)xmalloc(num_procs*sizeof(int)); + MPI_Gather(&local_lines,1,MPI_INT,all_lines,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(all_lines,num_procs,MPI_INT,0,MPI_COMM_WORLD); + all_cols=(int*)xmalloc(num_procs*sizeof(int)); + MPI_Gather(&local_cols,1,MPI_INT,all_cols,1,MPI_INT,0,MPI_COMM_WORLD); + MPI_Bcast(all_cols,num_procs,MPI_INT,0,MPI_COMM_WORLD); + + for(i=0;icols) cols=all_cols[i]; + for(i=0;iEnum())){ + element=(Element*)(*object); + //elemcols=element->InputToResult(&solution,enum_type); + } + } + + /*Allocate output pointer*/ + *psolution=solution; + + /* Free ressources: */ + xfree((void**)&all_cols); + xfree((void**)&all_lines); } /*}}}*/ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 3868) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 3869) @@ -93,4 +93,6 @@ void FieldAverageOntoVertices(Vec fieldsum,Vec connectivity,double* field); void FieldDepthAverageAtBase(Vec field,double* field_serial,char* fieldname); + void DepthAverageInputAtBase(int enum_type); + void InputExtrude(int enum_type); int DeleteObject(Object* object); void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type); @@ -99,4 +101,5 @@ int FindResult(void* pvalue, char* name); void FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse); + void InputToResult(Mat* psolution,int enum_type); void UpdateVertexPositions(double* thickness,double* bed); void OutputRifts(Vec riftproperties); @@ -128,4 +131,5 @@ /*numerics: {{{1*/ int AddInput(Input* in_input); + Input* GetInput(int enum_name); Inputs* SpawnTriaInputs(int* indices); Index: /issm/trunk/src/c/DataSet/Inputs.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3868) +++ /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3869) @@ -529,4 +529,21 @@ } /*}}}*/ +/*FUNCTION Inputs::GetInput{{{1*/ +Input* Inputs::GetInput(int enum_name){ + + vector::iterator object; + Input* input=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + input=(Input*)(*object); + + if (input->EnumType()==enum_name){ + return input; + } + } + ISSMERROR("input with enum %i (%s) not found",enum_name,EnumAsString(enum_name)); +} +/*}}}*/ /*FUNCTION Inputs::ChangeEnum{{{1*/ void Inputs::ChangeEnum(int oldenumtype,int newenumtype){ Index: /issm/trunk/src/c/DepthAverageInputx/DepthAverageInputx.cpp =================================================================== --- /issm/trunk/src/c/DepthAverageInputx/DepthAverageInputx.cpp (revision 3868) +++ /issm/trunk/src/c/DepthAverageInputx/DepthAverageInputx.cpp (revision 3869) @@ -11,13 +11,13 @@ void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type){ + + /*First, get elements*/ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); - - /*First, get elements and loads configured: */ - elements->Configure(elements,loads, nodes,vertices, materials,parameters); - nodes->Configure(elements,loads, nodes,vertices, materials,parameters); - loads->Configure(elements, loads, nodes,vertices, materials,parameters); - parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + /*First depth-average inputs at base of the glacier*/ + elements->DepthAverageInputAtBase(enum_type); - + /*Then extrude vertically the new inputs*/ + elements->InputExtrude(enum_type); } Index: /issm/trunk/src/c/ExtrudeInputx/ExtrudeInputx.cpp =================================================================== --- /issm/trunk/src/c/ExtrudeInputx/ExtrudeInputx.cpp (revision 3868) +++ /issm/trunk/src/c/ExtrudeInputx/ExtrudeInputx.cpp (revision 3869) @@ -12,5 +12,9 @@ void ExtrudeInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type){ - ISSMERROR(" not supported yet!"); + /*First, get elements*/ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Then extrude vertically the new inputs*/ + elements->InputExtrude(enum_type); } Index: /issm/trunk/src/c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.cpp =================================================================== --- /issm/trunk/src/c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.cpp (revision 3868) +++ /issm/trunk/src/c/FieldAverageOntoVerticesx/FieldAverageOntoVerticesx.cpp (revision 3869) @@ -23,5 +23,5 @@ Vec field=NULL; - /*Initioalize intermediary*/ + /*Initialize intermediary*/ found=parameters->FindParam(&numberofvertices,NumberOfVerticesEnum); connectivity=NewVec(numberofvertices); Index: /issm/trunk/src/c/InputToResultx/InputToResultx.cpp =================================================================== --- /issm/trunk/src/c/InputToResultx/InputToResultx.cpp (revision 3868) +++ /issm/trunk/src/c/InputToResultx/InputToResultx.cpp (revision 3869) @@ -13,5 +13,21 @@ void InputToResultx(Result** presult,DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,Parameters* parameters,int enum_type,int id, double time, int step){ - ISSMERROR(" not supported yet!"); + /*Intermediary output vector*/ + Mat solution=NULL; + Result* result=NULL; + + /*First, get elements*/ + elements->Configure(elements,loads, nodes,vertices, materials,parameters); + + /*Then extrude vertically the new inputs*/ + elements->InputToResult(&solution,enum_type); + + /*Create Result and clean up*/ + /*Mat support in result to be added!!!!!!*/ + //result=new Result(id,time,step,enum_type,solution); + MatFree(&solution); + + /*Assign output pointer*/ + *presult=result; } Index: /issm/trunk/src/c/objects/Elements/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Beam.h (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Beam.h (revision 3869) @@ -52,14 +52,14 @@ Object* copy(); void SetClone(int* minranks); - void UpdateInputsFromVector(double* vector, int name, int type); - void UpdateInputsFromVector(int* vector, int name, int type); - void UpdateInputsFromVector(bool* vector, int name, int type); - void UpdateInputsFromConstant(double constant, int name){ISSMERROR("Not implemented yet!");} - void UpdateInputsFromConstant(int constant, int name){ISSMERROR("Not implemented yet!");} - void UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} + void UpdateInputsFromVector(double* vector, int name, int type); + void UpdateInputsFromVector(int* vector, int name, int type); + void UpdateInputsFromVector(bool* vector, int name, int type); + void UpdateInputsFromConstant(double constant, int name){ISSMERROR("Not implemented yet!");} + void UpdateInputsFromConstant(int constant, int name){ISSMERROR("Not implemented yet!");} + void UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");} - void UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type); + void UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type); - + void DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");}; /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 3869) @@ -42,4 +42,5 @@ virtual double CostFunction(int analysis_type,int sub_analysis_type)=0; virtual double SurfaceArea(int analysis_type,int sub_analysis_type)=0; + virtual void DepthAverageInputAtBase(int enum_type)=0; virtual void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type)=0; virtual void ComputePressure(Vec p_g, int analysis_type,int sub_analysis_type)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 3869) @@ -337,4 +337,22 @@ } else return false; +} +/*}}}*/ +/*FUNCTION Penta::IsOnSurface{{{1*/ +bool Penta::IsOnSurface(void){ + + bool onsurface; + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + return onsurface; + +} +/*}}}*/ +/*FUNCTION Penta::GetUpperElement{{{1*/ +Penta* Penta::GetUpperElement(void){ + + Penta* upper_penta=NULL; + upper_penta=(Penta*)neighbors[1]; //first one under, second one above + return upper_penta; + } /*}}}*/ @@ -2962,4 +2980,46 @@ } //if (extrude) +} +/*}}}*/ +/*FUNCTION Penta::InputExtrude {{{1*/ +void Penta::InputExtrude(int enum_type){ + + bool onbed; + bool onsurface; + + Penta* penta=NULL; + + /*recover parameters: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + + /*Are we on the base, not on the surface?:*/ + if(onbed==1 & onsurface==0){ + + /*this element is on the bed. We are going to, follow the upper element until we reach the surface. At each upper element, + * we'll add the input of the lower element*/ + + penta=this; + + for(;;){ + + /* get upper Penta*/ + penta=penta->GetUpperElement(); + + /*Add input of the basal element*/ + penta->inputs->AddInput(this->inputs->GetInput(enum_type)); + + /*Stop if we have reached the surface*/ + if (penta->IsOnSurface()) break; + + } + } + + return; +} +/*}}}*/ +/*FUNCTION Penta::DepthAverageInputAtBase{{{1*/ +void Penta::DepthAverageInputAtBase(int enum_type){ + ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)"); } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 3869) @@ -50,5 +50,7 @@ int Enum(); int Id(); - bool IsInput(int name); + bool IsInput(int name); + bool IsOnSurface(void); + Penta* GetUpperElement(void); void Marshall(char** pmarshalled_dataset); int MarshallSize(); @@ -98,4 +100,6 @@ void GetNodalFunctions(double* l1l6, double* gauss_coord); void FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed); + void InputExtrude(int enum_type); + void DepthAverageInputAtBase(int enum_type); void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type); void ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/Elements/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Sing.h (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Sing.h (revision 3869) @@ -61,4 +61,5 @@ void UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type); + void DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");}; /*}}}*/ /*numerics: {{{1*/ Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3869) @@ -470,5 +470,27 @@ /*FUNCTION Tria::UpdateInputsFromSolutionPrognostic {{{1*/ void Tria::UpdateInputsFromSolutionPrognostic(double* solution, int analysis_type, int sub_analysis_type){ - ISSMERROR(" not supported yet!"); + + int i; + + const int numvertices=3; + const int numdofpervertex=1; + const int numdof=numdofpervertex*numvertices; + + int doflist[numdof]; + double values[numdof]; + double thickness[numvertices]; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iinputs->AddInput(new TriaVertexInput(ThicknessEnum,values)); } /*}}}*/ @@ -4663,2 +4685,33 @@ } /*}}}*/ +/*FUNCTION Tria::DepthAverageInputAtBase {{{1*/ +void Tria::DepthAverageInputAtBase(int enum_type){ + + /*New input*/ + Input* oldinput=NULL; + Input* newinput=NULL; + + /*copy input of enum_type*/ + oldinput=this->inputs->GetInput(enum_type); + newinput=(Input*)oldinput->copy(); + + /*Assign new name (average)*/ + switch(enum_type){ + + case VxEnum: + newinput->ChangeEnum(VxAverageEnum); + break; + + case VyEnum: + newinput->ChangeEnum(VyAverageEnum); + break; + + default: + ISSMERROR("enum_type: %i (%s) not suppoerted yet",enum_type,EnumAsString(enum_type)); + } + + /*Add new input to current element*/ + this->inputs->AddInput(newinput); + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 3868) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 3869) @@ -52,5 +52,5 @@ int MyRank(); void SetClone(int* minranks); - + void DepthAverageInputAtBase(int enum_type); /*}}}*/ /*FUNCTION element numerical routines {{{1*/