Changeset 4460

Timestamp:

07/08/10 11:02:13 (15 years ago)

Author:

Mathieu Morlighem

Message:

Began DepthAverageAtBase, to be continued

File:

• issm/trunk/src/c/objects/Elements/Penta.cpp (modified) (1 diff)

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

@@ -1074 +1074 @@
 /*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
 void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
-        ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)");
+        ISSMERROR("Not implemented yet (see Node::FieldDepthAverageAtBase)");
+
+        /*Intermediaries*/
+        const int numvertices=6;
+        bool onbed;
+        bool onsurface;
+        int  i;
+
+        Penta* penta=NULL;
+        Input* original_input=NULL;
+        Input* integrated_input=NULL;
+        Input* element_thickness_input=NULL;
+        Input* total_thickness_input=NULL;
+        Input* depth_averaged_input=NULL;
+
+        double  xyz_list[numvertices][3];
+        double  Helem_list[numvertices];
+        double  zeros_list[numvertices]={0.0};
+
+        /*recover parameters: */
+        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+        /*Are we on the base, not on the surface?:*/
+        if(onbed==1){
+
+                /*OK, we are on bed. Initialize inputs*/
+                total_thickness_input=new PentaVertexInput(ThicknessEnum,zeros_list);
+                depth_averaged_input =new PentaVertexInput(average_enum_type,zeros_list);
+
+                /*Now follow all the upper element from the base to the surface to integrate the input*/
+                penta=this;
+                for(;;){
+
+                        /*Step1: Get original input (to be depth avegaged): */
+                        original_input=(Input*)penta->inputs->GetInput(enum_type);
+                        if(!original_input) ISSMERROR("%s%s"," could not find input with enum:",EnumAsString(enum_type));
+
+                        /*Step2: Create element thickness input*/
+                        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+                        for(i=0;i<3;i++){
+                                Helem_list[i]=xyz_list[i+3][2]-xyz_list[i][2];
+                                Helem_list[i+3]=Helem_list[i];
+                        }
+                        element_thickness_input=new PentaVertexInput(ThicknessEnum,Helem_list);
+
+                        /*Step3: Vertically integrate original input and update totalthickness_input*/
+                        //integrated_input=original_input->VerticallyIntegrate(element_thickness_input); //TO BE ADDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+                        total_thickness_input->AXPY(element_thickness_input,1.0);
+
+                        /*Stop if we have reached the surface, else, take upper penta*/
+                        if (penta->IsOnSurface()){
+                                break;
+                        }
+                        else{
+                                /* get upper Penta*/
+                                penta=penta->GetUpperElement();
+                                ISSMASSERT(penta->Id()!=this->id);
+                        }
+                }
+
+                /*OK, now we only need to divide the depth integrated input by the total thickness!*/
+                //depth_averaged_input=integrated_input->InputPointwiseDivide(total_thickness_input);  //TO BE ADDED !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+                depth_averaged_input->ChangeEnum(average_enum_type);
+
+                /*Finally, add to inputs*/
+                this->inputs->AddInput((Input*)depth_averaged_input);
+        }
+
+        return;
 }
 /*}}}*/