Changeset 4931


Ignore:
Timestamp:
08/02/10 14:55:29 (15 years ago)
Author:
Mathieu Morlighem
Message:

Now Matice uses B2d in 2d and B in 3d. B2d is added to matice's inputs when penta spawns a tria

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

Legend:

Unmodified
Added
Removed

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

    r4573 r4931  
    1919        for (i=0;i<elements->Size();i++){
    2020                element=(Element*)elements->GetObjectByOffset(i);
    21                 element->InputDepthAverageAtBase(enum_type,average_enum_type);
     21                element->InputDepthAverageAtBase(enum_type,average_enum_type,ElementsEnum);
    2222        }
    2323

  • issm/trunk/src/c/objects/Elements/Element.h

    r4880 r4931  
    4343                virtual double CostFunction(void)=0;
    4444                virtual double SurfaceArea(void)=0;
    45                 virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type)=0;
     45                virtual void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum)=0;
    4646                virtual void   ComputeBasalStress(Vec sigma_b)=0;
    4747                virtual void   ComputePressure(Vec p_g)=0;

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

    r4929 r4931  
    1414#include <string.h>
    1515#include "../objects.h"
    16 #include "../../EnumDefinitions/EnumDefinitions.h"
    1716#include "../../shared/shared.h"
    1817#include "../../include/include.h"
    … …  
    10281027/*}}}*/
    10291028/*FUNCTION Penta::InputDepthAverageAtBase{{{1*/
    1030 void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){
     1029void  Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
    10311030
    10321031        /*Intermediaries*/
    10331032        const int numvertices=6;
    10341033        bool onbed;
    1035         int  i;
     1034        int  step,i;
    10361035
    10371036        Penta* penta=NULL;
    … …  
    10551054        /*OK, we are on bed. Initialize global inputs as 0*/
    10561055        total_thickness_input =new PentaVertexInput(ThicknessEnum,zeros_list);
    1057         total_integrated_input=new PentaVertexInput(average_enum_type,zeros_list);
    10581056
    10591057        /*Now follow all the upper element from the base to the surface to integrate the input*/
    10601058        penta=this;
     1059        step =0;
    10611060        for(;;){
    10621061
    10631062                /*Step1: Get original input (to be depth avegaged): */
    1064                 original_input=(Input*)penta->inputs->GetInput(enum_type);
     1063                if (object_enum==ElementsEnum)
     1064                 original_input=(Input*)penta->inputs->GetInput(enum_type);
     1065                else if (object_enum==MaterialsEnum)
     1066                 original_input=(Input*)penta->matice->inputs->GetInput(enum_type);
     1067                else
     1068                 ISSMERROR("object %s not supported yet",EnumAsString(object_enum));
    10651069                if(!original_input) ISSMERROR("could not find input with enum %s",EnumAsString(enum_type));
     1070
     1071                /*If first time, initialize total_integrated_input*/
     1072                if (step==0){
     1073                        if (original_input->Enum()==PentaVertexInputEnum)
     1074                         total_integrated_input=new PentaVertexInput(average_enum_type,zeros_list);
     1075                        else if (original_input->Enum()==DoubleInputEnum)
     1076                         total_integrated_input=new DoubleInput(average_enum_type,0.0);
     1077                        else
     1078                         ISSMERROR("object %s not supported yet",EnumAsString(original_input->Enum()));
     1079                }
    10661080
    10671081                /*Step2: Create element thickness input*/
    … …  
    10861100
    10871101                /*Stop if we have reached the surface, else, take upper penta*/
    1088                 if (penta->IsOnSurface()){
    1089                         break;
    1090                 }
    1091                 else{
    1092                         /* get upper Penta*/
    1093                         penta=penta->GetUpperElement();
    1094                         ISSMASSERT(penta->Id()!=this->id);
    1095                 }
     1102                if (penta->IsOnSurface()) break;
     1103
     1104                /* get upper Penta*/
     1105                penta=penta->GetUpperElement();
     1106                ISSMASSERT(penta->Id()!=this->id);
     1107
     1108                /*increase couter*/
     1109                step++;
    10961110        }
    10971111
    … …  
    11051119
    11061120        /*Finally, add to inputs*/
    1107         this->inputs->AddInput((Input*)depth_averaged_input);
    1108         return;
     1121        if (object_enum==ElementsEnum)
     1122         this->inputs->AddInput((Input*)depth_averaged_input);
     1123        else if (object_enum==MaterialsEnum)
     1124         this->matice->inputs->AddInput((Input*)depth_averaged_input);
     1125        else
     1126         ISSMERROR("object %s not supported yet",EnumAsString(object_enum));
    11091127}
    11101128/*}}}*/
    … …  
    21302148                /*This element should be collapsed into a tria element at its base. Create this tria element,
    21312149                 *and use its CreateKMatrix functionality to fill the global stiffness matrix: */
     2150
     2151                /*Depth Averaging B*/
     2152                this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
     2153
     2154                /*Call Tria function*/
    21322155                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    21332156                tria->CreateKMatrix(Kgg);
    21342157                delete tria;
     2158
     2159                /*Delete B averaged*/
     2160                this->matice->inputs->DeleteInput(RheologyB2dEnum);
     2161
    21352162                return;
    21362163        }
    … …  
    39413968
    39423969        /*Delete Vx and Vy averaged*/
    3943         this-inputs->DeleteInput(VxAverageEnum);
    3944         this-inputs->DeleteInput(VyAverageEnum);
     3970        this->inputs->DeleteInput(VxAverageEnum);
     3971        this->inputs->DeleteInput(VyAverageEnum);
    39453972
    39463973        return;
    … …  
    52285255        this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,values));
    52295256
    5230         /*Also update the rheology WARNING: ONLY FOR PATTYN and STOKES FOR NOW*/
     5257        /*Also update the rheology*/
    52315258        inputs->GetParameterValue(&collapse,CollapseEnum);
    5232         if (!collapse) this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
     5259        this->matice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B));
    52335260
    52345261}

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

    r4927 r4931  
    2222#include "../../include/include.h"
    2323#include "../../shared/Exceptions/exceptions.h"
     24#include "../../EnumDefinitions/EnumDefinitions.h"
    2425/*}}}*/
    2526
    … …  
    8889                void   InputControlConstrain(int control_type,double cm_min, double cm_max);
    8990                bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    90                 void   InputDepthAverageAtBase(int enum_type,int average_enum_type);
     91                void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=ElementsEnum);
    9192                void   InputDuplicate(int original_enum,int new_enum);
    9293                void   InputScale(int enum_type,double scale_factor);

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

    r4927 r4931  
    1414#include <string.h>
    1515#include "../objects.h"
    16 #include "../../EnumDefinitions/EnumDefinitions.h"
    1716#include "../../shared/shared.h"
    1817#include "../../Container/Container.h"
    … …  
    629628                }
    630629                else if (control_type==RheologyBEnum){
    631                         inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBEnum);
     630                        inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyB2dEnum);
    632631                        Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight;
    633632                }
    … …  
    870869                GradjDrag(gradient);
    871870        }
    872         else if (control_type==RheologyBEnum){
     871        else if (control_type==RheologyB2dEnum){
    873872                GradjB(gradient);
    874873        }
    … …  
    966965        adjointx_input=inputs->GetInput(AdjointxEnum);
    967966        adjointy_input=inputs->GetInput(AdjointyEnum);
    968         rheologyb_input=inputs->GetInput(RheologyBEnum);
     967        rheologyb_input=inputs->GetInput(RheologyB2dEnum);
    969968
    970969        /* Start  looping on the number of gaussian points: */
    … …  
    12311230
    12321231        /*Move input to Material if required (needed if control method) TO BE IMPROVED*/
    1233         if (YEnum==RheologyBEnum){
     1232        if (YEnum==RheologyBEnum || YEnum==RheologyB2dEnum){
    12341233                this->matice->inputs->AddInput((Input*)yinput->copy());
    12351234        }
    … …  
    12871286/*}}}*/
    12881287/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/
    1289 void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){
     1288void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){
    12901289
    12911290        /*New input*/
    … …  
    12941293
    12951294        /*copy input of enum_type*/
    1296         oldinput=this->inputs->GetInput(enum_type);
     1295        if (object_enum==ElementsEnum)
     1296         oldinput=(Input*)this->inputs->GetInput(enum_type);
     1297        else if (object_enum==MaterialsEnum)
     1298         oldinput=(Input*)this->matice->inputs->GetInput(enum_type);
     1299        else
     1300         ISSMERROR("object %s not supported yet",EnumAsString(object_enum));
    12971301        if(!oldinput)ISSMERROR("%s%s"," could not find old input with enum: ",EnumAsString(enum_type));
    12981302        newinput=(Input*)oldinput->copy();
    … …  
    13021306
    13031307        /*Add new input to current element*/
    1304         this->inputs->AddInput(newinput);
     1308        if (object_enum==ElementsEnum)
     1309         this->inputs->AddInput((Input*)newinput);
     1310        else if (object_enum==MaterialsEnum)
     1311         this->matice->inputs->AddInput((Input*)newinput);
     1312        else
     1313         ISSMERROR("object %s not supported yet",EnumAsString(object_enum));
    13051314
    13061315}
    … …  
    21262135        if (iomodel->rheology_B) {
    21272136                for(i=0;i<3;i++)nodeinputs[i]=iomodel->rheology_B[tria_vertex_ids[i]-1];
    2128                 this->inputs->AddInput(new TriaVertexInput(RheologyBEnum,nodeinputs));
     2137                this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs));
    21292138        }
    21302139        if (iomodel->control_parameter) {
    … …  
    60486057                                name==VyEnum ||
    60496058                                name==RheologyBEnum ||
     6059                                name==RheologyB2dEnum ||
    60506060                                name==RheologyNEnum ||
    60516061                                name==FitEnum ||

  • issm/trunk/src/c/objects/Elements/Tria.h

    r4927 r4931  
    2020#include "../../include/include.h"
    2121#include "../../shared/Exceptions/exceptions.h"
     22#include "../../EnumDefinitions/EnumDefinitions.h"
    2223/*}}}*/
    2324
    … …  
    8485                void   InputControlConstrain(int control_type,double cm_min, double cm_max);
    8586                bool   InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
    86                 void   InputDepthAverageAtBase(int enum_type,int average_enum_type);
     87                void   InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum=ElementsEnum);
    8788                void   InputDuplicate(int original_enum,int new_enum);
    8889                void   InputScale(int enum_type,double scale_factor);

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

    r4927 r4931  
    2323}
    2424/*}}}*/
    25 /*FUNCTION DoubleInput::DoubleInput(double* values){{{1*/
     25/*FUNCTION DoubleInput::DoubleInput(double value){{{1*/
    2626DoubleInput::DoubleInput(int in_enum_type,IssmDouble in_value){
    2727
    … …  
    253253}
    254254/*}}}*/
     255/*FUNCTION DoubleInput::VerticallyIntegrate{{{1*/
     256void DoubleInput::VerticallyIntegrate(Input* thickness_input){
     257
     258        /*Intermediaries*/
     259        double thickness_value;
     260
     261        /*Check that input provided is a thickness*/
     262        if (thickness_input->EnumType()!=ThicknessEnum) ISSMERROR("Input provided is not a Thickness (enum_type is %s)",EnumAsString(thickness_input->EnumType()));
     263
     264        /*vertically integrate depending on type:*/
     265        switch(thickness_input->Enum()){
     266
     267                case PentaVertexInputEnum:
     268                        thickness_input->GetParameterAverage(&thickness_value);
     269                        this->value=this->value*thickness_value;
     270                        return;
     271
     272                default:
     273                        ISSMERROR("not implemented yet");
     274        }
     275}
     276/*}}}*/
     277/*FUNCTION DoubleInput::PointwiseDivide{{{1*/
     278Input* DoubleInput::PointwiseDivide(Input* inputB){
     279
     280        /*Ouput*/
     281        DoubleInput* outinput=NULL;
     282
     283        /*Intermediaries*/
     284        double       AdotBvalue;
     285        double       Bvalue;
     286
     287        /*Check that inputB is of the same type*/
     288        inputB->GetParameterAverage(&Bvalue);
     289
     290        /*Create new DoubleInput*/
     291        outinput=new DoubleInput(this->enum_type,this->value/Bvalue);
     292
     293        /*Return output pointer*/
     294        return outinput;
     295
     296}
     297/*}}}*/

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

    r4927 r4931  
    3838                int   EnumType();
    3939                Input* SpawnTriaInput(int* indices);
    40                 Input* PointwiseDivide(Input* inputB){ISSMERROR("not implemented yet");};
     40                Input* PointwiseDivide(Input* inputB);
    4141                ElementResult* SpawnResult(int step, double time);
    4242                /*}}}*/
    … …  
    6262                void Constrain(double cm_min, double cm_max);
    6363                void Extrude(void){ISSMERROR("not supported yet");};
    64                 void VerticallyIntegrate(Input* thickness_input){ISSMERROR("not supported yet");};
     64                void VerticallyIntegrate(Input* thickness_input);
    6565                void GetVectorFromInputs(Vec vector,int* doflist);
    6666                void GetValuesPtr(double** pvalues,int* pnum_values);

  • issm/trunk/src/c/objects/Materials/Matice.cpp

    r4920 r4931  
    4747                if (iomodel->rheology_B) {
    4848                        for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->rheology_B[int(iomodel->elements[num_vertices*index+i]-1)];
    49                         this->inputs->AddInput(new TriaVertexInput(RheologyBEnum,nodeinputs));
     49                        this->inputs->AddInput(new TriaVertexInput(RheologyB2dEnum,nodeinputs));
    5050                }
    5151
    … …  
    245245}
    246246/*}}}*/
     247/*FUNCTION Matice::GetB2d {{{1*/
     248double Matice::GetB2d(){
     249
     250        /*Output*/
     251        double B2d;
     252
     253        inputs->GetParameterAverage(&B2d,RheologyB2dEnum);
     254        return B2d;
     255}
     256/*}}}*/
    247257/*FUNCTION Matice::GetN {{{1*/
    248258double Matice::GetN(){
    … …  
    280290
    281291        /*Get B and n*/
    282         B=GetB();
     292        B=GetB2d();
    283293        n=GetN();
    284294
    … …  
    488498
    489499        /*Get B and n*/
    490         B=GetB();
     500        B=GetB2d();
    491501        n=GetN();
    492502

  • issm/trunk/src/c/objects/Materials/Matice.h

    r4920 r4931  
    5858                void   GetViscosityComplement(double* pviscosity_complement, double* pepsilon);
    5959                double GetB();
     60                double GetB2d();
    6061                double GetN();
    6162                /*}}}*/
Note: See TracChangeset for help on using the changeset viewer.