Changeset 4471


Ignore:
Timestamp:
07/08/10 14:58:01 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added DepthAverageInput module

Location:
issm/trunk/src/c
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/modules/InputDepthAveragex/InputDepthAveragex.cpp

    r4218 r4471  
    2525        }
    2626
    27 
    2827        /*Then extrude vertically the new inputs*/
    29         InputExtrudex( elements,nodes,vertices,loads,materials,parameters,enum_type);
     28        InputExtrudex( elements,nodes,vertices,loads,materials,parameters,average_enum_type);
    3029}

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4470 r4471  
    10801080/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
    10811081void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
    1082         ISSMERROR("Not implemented yet (see Node::FieldDepthAverageAtBase)");
    10831082
    10841083        /*Intermediaries*/
    … …  
    10901089        Penta* penta=NULL;
    10911090        Input* original_input=NULL;
    1092         Input* integrated_input=NULL;
     1091        Input* element_integrated_input=NULL;
     1092        Input* total_integrated_input=NULL;
    10931093        Input* element_thickness_input=NULL;
    10941094        Input* total_thickness_input=NULL;
    … …  
    11021102        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
    11031103
    1104         /*Are we on the base, not on the surface?:*/
    1105         if(onbed==1){
    1106 
    1107                 /*OK, we are on bed. Initialize inputs*/
    1108                 total_thickness_input=new PentaVertexInput(ThicknessEnum,zeros_list);
    1109                 depth_averaged_input =new PentaVertexInput(average_enum_type,zeros_list);
    1110 
    1111                 /*Now follow all the upper element from the base to the surface to integrate the input*/
    1112                 penta=this;
    1113                 for(;;){
    1114 
    1115                         /*Step1: Get original input (to be depth avegaged): */
    1116                         original_input=(Input*)penta->inputs->GetInput(enum_type);
    1117                         if(!original_input) ISSMERROR("%s%s"," could not find input with enum:",EnumAsString(enum_type));
    1118 
    1119                         /*Step2: Create element thickness input*/
    1120                         GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
    1121                         for(i=0;i<3;i++){
    1122                                 Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2];
    1123                                 Helem_list[i+3]=Helem_list[i];
    1124                         }
    1125                         element_thickness_input=new PentaVertexInput(ThicknessEnum,Helem_list);
    1126 
    1127                         /*Step3: Vertically integrate original input and update totalthickness_input*/
    1128                         //integrated_input=original_input->VerticallyIntegrate(element_thickness_input); //TO BE ADDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1129                         total_thickness_input->AXPY(element_thickness_input,1.0);
    1130 
    1131                         /*Stop if we have reached the surface, else, take upper penta*/
    1132                         if (penta->IsOnSurface()){
    1133                                 break;
    1134                         }
    1135                         else{
    1136                                 /* get upper Penta*/
    1137                                 penta=penta->GetUpperElement();
    1138                                 ISSMASSERT(penta->Id()!=this->id);
    1139                         }
     1104        /*Are we on the base? If not, return*/
     1105        if(!onbed) return;
     1106
     1107        /*OK, we are on bed. Initialize global inputs as 0*/
     1108        total_thickness_input =new PentaVertexInput(ThicknessEnum,zeros_list);
     1109        total_integrated_input=new PentaVertexInput(average_enum_type,zeros_list);
     1110
     1111        /*Now follow all the upper element from the base to the surface to integrate the input*/
     1112        penta=this;
     1113        for(;;){
     1114
     1115                /*Step1: Get original input (to be depth avegaged): */
     1116                original_input=(Input*)penta->inputs->GetInput(enum_type);
     1117                if(!original_input) ISSMERROR("could not find input with enum %s",EnumAsString(enum_type));
     1118
     1119                /*Step2: Create element thickness input*/
     1120                GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
     1121                for(i=0;i<3;i++){
     1122                        Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2];
     1123                        Helem_list[i+3]=Helem_list[i];
    11401124                }
    1141 
    1142                 /*OK, now we only need to divide the depth integrated input by the total thickness!*/
    1143                 //depth_averaged_input=integrated_input->InputPointwiseDivide(total_thickness_input);  //TO BE ADDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1144                 depth_averaged_input->ChangeEnum(average_enum_type);
    1145 
    1146                 /*Finally, add to inputs*/
    1147                 this->inputs->AddInput((Input*)depth_averaged_input);
    1148         }
    1149 
     1125                element_thickness_input=new PentaVertexInput(ThicknessEnum,Helem_list);
     1126
     1127                /*Step3: Vertically integrate A COPY of the original*/
     1128                element_integrated_input=(Input*)original_input->copy();
     1129                element_integrated_input->VerticallyIntegrate(element_thickness_input);
     1130
     1131                /*Add contributions to global inputs*/
     1132                total_integrated_input->AXPY(element_integrated_input,1.0);
     1133                total_thickness_input ->AXPY(element_thickness_input,1.0);
     1134
     1135                /*Clean up*/
     1136                delete element_thickness_input;
     1137                delete element_integrated_input;
     1138
     1139                /*Stop if we have reached the surface, else, take upper penta*/
     1140                if (penta->IsOnSurface()){
     1141                        break;
     1142                }
     1143                else{
     1144                        /* get upper Penta*/
     1145                        penta=penta->GetUpperElement();
     1146                        ISSMASSERT(penta->Id()!=this->id);
     1147                }
     1148        }
     1149
     1150        /*OK, now we only need to divide the depth integrated input by the total thickness!*/
     1151        depth_averaged_input=total_integrated_input->PointwiseDivide(total_thickness_input);
     1152        depth_averaged_input->ChangeEnum(average_enum_type);
     1153
     1154        /*Clean up*/
     1155        delete total_thickness_input;
     1156        delete total_integrated_input;
     1157
     1158        /*Finally, add to inputs*/
     1159        this->inputs->AddInput((Input*)depth_averaged_input);
    11501160        return;
    11511161}

  • issm/trunk/src/c/objects/Inputs/BeamVertexInput.h

    r4274 r4471  
    4141                Input* SpawnBeamInput(int* indices);
    4242                Input* SpawnTriaInput(int* indices);
     43                Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
    4344                ElementResult* SpawnResult(int step, double time);
    4445                /*}}}*/
    … …  
    7071                void Constrain(double cm_min, double cm_max);
    7172                void Extrude(void){ISSMERROR("not supported yet");};
     73                void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
    7274                void GetVectorFromInputs(Vec vector,int* doflist);
    7375                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Inputs/BoolInput.h

    r4274 r4471  
    4141                Input* SpawnBeamInput(int* indices);
    4242                Input* SpawnTriaInput(int* indices);
     43                Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
    4344                ElementResult* SpawnResult(int step, double time);
    4445                /*}}}*/
    … …  
    7071                void Constrain(double cm_min, double cm_max);
    7172                void Extrude(void){ISSMERROR("not supported yet");};
     73                void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
    7274                void GetVectorFromInputs(Vec vector,int* doflist);
    7375                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Inputs/DoubleInput.h

    r4274 r4471  
    4040                Input* SpawnBeamInput(int* indices);
    4141                Input* SpawnTriaInput(int* indices);
     42                Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
    4243                ElementResult* SpawnResult(int step, double time);
    4344                /*}}}*/
    … …  
    6970                void Constrain(double cm_min, double cm_max);
    7071                void Extrude(void){ISSMERROR("not supported yet");};
     72                void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
    7173                void GetVectorFromInputs(Vec vector,int* doflist);
    7274                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Inputs/Input.h

    r4274 r4471  
    4747                virtual Input* SpawnBeamInput(int* indices)=0;
    4848                virtual Input* SpawnTriaInput(int* indices)=0;
     49                virtual Input* PointwiseDivide(Input* inputB)=0;
    4950                virtual ElementResult* SpawnResult(int step, double time)=0;
    5051                virtual void SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0;
    … …  
    5253                virtual void AXPY(Input* xinput,double scalar)=0;
    5354                virtual void Constrain(double cm_min, double cm_max)=0;
     55                virtual void VerticallyIntegrate(Input* thickness_input)=0;
    5456                virtual void Extrude()=0;
    5557                virtual void GetVectorFromInputs(Vec vector,int* doflist)=0;

  • issm/trunk/src/c/objects/Inputs/IntInput.h

    r4274 r4471  
    4141                Input* SpawnBeamInput(int* indices);
    4242                Input* SpawnTriaInput(int* indices);
     43                Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
    4344                ElementResult* SpawnResult(int step, double time);
    4445                /*}}}*/
    … …  
    7071                void Constrain(double cm_min, double cm_max);
    7172                void Extrude(void){ISSMERROR("not supported yet");};
     73                void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
    7274                void GetVectorFromInputs(Vec vector,int* doflist);
    7375                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp

    r4274 r4471  
    878878}
    879879/*}}}*/
    880 /*FUNCTION PentaVertexInput::SquareMin(double* psquaremin, bool process_units){{{1*/
     880/*FUNCTION PentaVertexInput::SquareMin{{{1*/
    881881void PentaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
    882882
    … …  
    901901}
    902902/*}}}*/
    903 /*FUNCTION PentaVertexInput::Scale(double scale_factor){{{1*/
     903/*FUNCTION PentaVertexInput::Scale{{{1*/
    904904void PentaVertexInput::Scale(double scale_factor){
    905905       
    … …  
    910910}
    911911/*}}}*/
    912 /*FUNCTION PentaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
     912/*FUNCTION PentaVertexInput::AXPY{{{1*/
    913913void PentaVertexInput::AXPY(Input* xinput,double scalar){
    914914
    … …  
    933933}
    934934/*}}}*/
    935 /*FUNCTION PentaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
     935/*FUNCTION PentaVertexInput::Constrain{{{1*/
    936936void PentaVertexInput::Constrain(double cm_min, double cm_max){
    937937
    … …  
    944944}
    945945/*}}}*/
    946 /*FUNCTION PentaVertexInput::Extrude(double cm_min, double cm_max){{{1*/
     946/*FUNCTION PentaVertexInput::Extrude{{{1*/
    947947void PentaVertexInput::Extrude(void){
    948948
    … …  
    954954}
    955955/*}}}*/
     956/*FUNCTION PentaVertexInput::VerticallyIntegrate{{{1*/
     957void PentaVertexInput::VerticallyIntegrate(Input* thickness_input){
     958
     959        /*Intermediaries*/
     960        int i;
     961        const int  numgrids = 6;
     962        int        num_thickness_values;
     963        double    *thickness_values = NULL;
     964
     965        /*Check that input provided is a thickness*/
     966        if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumAsString(thickness_input->EnumType()));
     967
     968        /*Get Thickness value pointer*/
     969        thickness_input->GetValuesPtr(&thickness_values,&num_thickness_values);
     970
     971        /*vertically integrate depending on type:*/
     972        switch(thickness_input->Enum()){
     973
     974                case PentaVertexInputEnum:
     975                        for(i=0;i<3;i++){
     976                                this->values[i]=0.5*(this->values[i]+this->values[i+3]) * thickness_values[i];
     977                                this->values[i+3]=this->values[i];
     978                        }
     979                        return;
     980
     981                default:
     982                        ISSMERROR("not implemented yet");
     983        }
     984}
     985/*}}}*/
     986/*FUNCTION PentaVertexInput::PointwiseDivide{{{1*/
     987Input* PentaVertexInput::PointwiseDivide(Input* inputB){
     988
     989        /*Ouput*/
     990        PentaVertexInput* outinput=NULL;
     991
     992        /*Intermediaries*/
     993        int               i;
     994        PentaVertexInput *xinputB     = NULL;
     995        int               B_numvalues;
     996        double           *B_values    = NULL;
     997        const int         numgrids    = 6;
     998        double            AdotBvalues[numgrids];
     999
     1000        /*Check that inputB is of the same type*/
     1001        if (inputB->Enum()!=PentaVertexInputEnum) ISSMERROR("Operation not permitted because inputB is of type %s",EnumAsString(inputB->Enum()));
     1002        xinputB=(PentaVertexInput*)inputB;
     1003
     1004        /*Create point wise sum*/
     1005        for(i=0;i<numgrids;i++){
     1006                ISSMASSERT(xinputB->values[i]!=0);
     1007                AdotBvalues[i]=this->values[i]/xinputB->values[i];
     1008        }
     1009
     1010        /*Create new Sing input (copy of current input)*/
     1011        outinput=new PentaVertexInput(this->enum_type,&AdotBvalues[0]);
     1012
     1013        /*Return output pointer*/
     1014        return outinput;
     1015
     1016}
     1017/*}}}*/
    9561018/*FUNCTION PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
    9571019void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){

  • issm/trunk/src/c/objects/Inputs/PentaVertexInput.h

    r4274 r4471  
    4040                Input* SpawnBeamInput(int* indices);
    4141                Input* SpawnTriaInput(int* indices);
     42                Input* PointwiseDivide(Input* inputB);
    4243                ElementResult* SpawnResult(int step, double time);
    4344                /*}}}*/
    … …  
    7980                void Constrain(double cm_min, double cm_max);
    8081                void Extrude(void);
     82                void VerticallyIntegrate(Input* thickness_input);
    8183                void GetVectorFromInputs(Vec vector,int* doflist);
    8284                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Inputs/SingVertexInput.h

    r4274 r4471  
    4040                Input* SpawnBeamInput(int* indices);
    4141                Input* SpawnTriaInput(int* indices);
     42                Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
    4243                ElementResult* SpawnResult(int step, double time);
    4344                /*}}}*/
    … …  
    6970                void Constrain(double cm_min, double cm_max);
    7071                void Extrude(void){ISSMERROR("not supported yet");};
     72                void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
    7173                void GetVectorFromInputs(Vec vector,int* doflist);
    7274                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Inputs/TriaVertexInput.h

    r4274 r4471  
    4141                Input* SpawnBeamInput(int* indices);
    4242                Input* SpawnTriaInput(int* indices);
     43                Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
    4344                ElementResult* SpawnResult(int step, double time);
    4445                /*}}}*/
    … …  
    7778                void Constrain(double cm_min, double cm_max);
    7879                void Extrude(void){ISSMERROR("not supported yet");};
     80                void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
    7981                void GetVectorFromInputs(Vec vector,int* doflist);
    8082                void GetValuesPtr(double** pvalues,int* pnum_values);
Note: See TracChangeset for help on using the changeset viewer.