Changeset 4832


Ignore:
Timestamp:
07/27/10 11:41:32 (15 years ago)
Author:
Mathieu Morlighem
Message:

Extrudion of Vx and Vy is now done by the elements

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

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4804 r4832  
    21252125/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
    21262126void  Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
     2127
     2128        /*Inputs*/
     2129        bool    collapsed,onbed;
     2130
     2131        /*Recover inputs*/
     2132        inputs->GetParameterValue(&collapsed,CollapseEnum);
     2133        inputs->GetParameterValue(&onbed,ElementOnBedEnum);
     2134
     2135        /*MacAyeal, everything is done by the element on bed*/
     2136        if (collapsed){
     2137                if (!onbed){
     2138                        return;
     2139                }
     2140                else{
     2141                        InputUpdateFromSolutionDiagnosticMacAyeal(solution);
     2142                        return;
     2143                }
     2144        }
     2145        else{
     2146                InputUpdateFromSolutionDiagnosticPattyn(solution);
     2147        }
     2148}
     2149
     2150/*}}}*/
     2151/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{1*/
     2152void  Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){
     2153
     2154        int i;
     2155
     2156        const int    numvertices=6;
     2157        const int    numdofpervertex=2;
     2158        const int    numdof=numdofpervertex*numvertices;
     2159
     2160        int          doflist[numdof];
     2161        double       values[numdof];
     2162        double       vx[numvertices];
     2163        double       vy[numvertices];
     2164        double       vz[numvertices];
     2165        double       vel[numvertices];
     2166        int          dummy;
     2167        double       pressure[numvertices];
     2168        double       surface[numvertices];
     2169        double       rho_ice,g;
     2170        double       xyz_list[numvertices][3];
     2171        double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
     2172
     2173        Input  *VzInput        = NULL;
     2174        double *VzPtr          = NULL;
     2175        Penta  *penta          = NULL;
     2176
     2177        /*OK, we are on bed. Now recover results*/
     2178
     2179        /*Get dof list: */
     2180        GetDofList(&doflist[0],&dummy);
     2181
     2182        /*Get node data: */
     2183        GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
     2184
     2185        /*Use the dof list to index into the solution vector: */
     2186        for(i=0;i<numdof;i++){
     2187                values[i]=solution[doflist[i]];
     2188        }
     2189
     2190        /*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
     2191        for(i=0;i<3;i++){
     2192                vx[i]  =values[i*numdofpervertex+0];
     2193                vy[i]  =values[i*numdofpervertex+1];
     2194                vx[i+3]=vx[i];
     2195                vy[i+3]=vy[i];
     2196        }
     2197
     2198        /*Get parameters fro pressure computation*/
     2199        rho_ice=matpar->GetRhoIce();
     2200        g=matpar->GetG();
     2201
     2202        /*Start looping over all elements above current element and update all inputs*/
     2203        penta=this;
     2204        for(;;){
     2205
     2206                /*Now Compute vel*/
     2207                VzInput=inputs->GetInput(VzEnum);
     2208                if (VzInput){
     2209                        if (VzInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumAsString(VzInput->Enum()));
     2210                        VzInput->GetValuesPtr(&VzPtr,&dummy);
     2211                        for(i=0;i<numvertices;i++) vz[i]=VzPtr[i];
     2212                }
     2213                else{for(i=0;i<numvertices;i++) vz[i]=0.0;}
     2214                for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
     2215
     2216                /*Now compute pressure*/
     2217                inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
     2218                for(i=0;i<numvertices;i++){
     2219                        pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     2220                }
     2221
     2222                /*Now, we have to move the previous Vx and Vy inputs  to old
     2223                 * status, otherwise, we'll wipe them off: */
     2224                penta->inputs->ChangeEnum(VxEnum,VxOldEnum);
     2225                penta->inputs->ChangeEnum(VyEnum,VyOldEnum);
     2226                penta->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
     2227
     2228                /*Add vx and vy as inputs to the tria element: */
     2229                penta->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
     2230                penta->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
     2231                penta->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
     2232                penta->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
     2233
     2234                /*Stop if we have reached the surface*/
     2235                if (penta->IsOnSurface()) break;
     2236
     2237                /* get upper Penta*/
     2238                penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id);
     2239        }
     2240}
     2241/*}}}*/
     2242/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{1*/
     2243void  Penta::InputUpdateFromSolutionDiagnosticPattyn(double* solution){
    21272244       
    21282245        int i;
     
    21452262        double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
    21462263       
    2147         Input*       VzInput=NULL;
    2148         double*      VzPtr=NULL;
     2264        Input  *VzInput        = NULL;
     2265        double *VzPtr          = NULL;
    21492266
    21502267        /*Get dof list: */
     
    21862303        g=matpar->GetG();
    21872304        inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
    2188        
     2305
    21892306        for(i=0;i<numvertices;i++){
    21902307                pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
     
    21992316        this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
    22002317        this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
    2201         this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
     2318        this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
    22022319        this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
    22032320}
     
    49165033
    49175034        /*Are we on the base, not on the surface?:*/
    4918         if(onbed==1){
     5035        if(onbed){
    49195036
    49205037                /*OK, we are on bed. we will follow the steps:
  • issm/trunk/src/c/objects/Elements/Penta.h

    r4780 r4832  
    177177                void    InputUpdateFromSolutionBalancedvelocities( double* solutiong);
    178178                void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
     179                void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
     180                void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
    179181                void    InputUpdateFromSolutionDiagnosticHutter( double* solutiong);
    180182                void    InputUpdateFromSolutionDiagnosticVert( double* solutiong);
  • issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp

    r4592 r4832  
    128128        }
    129129       
    130 
    131         /*extrude if we are in 3D: */
    132         if (dim==3){
    133                 if(verbose)_printf_("%s\n","extruding velocity of collapsed elements...");
    134                 InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,true); //true means that the velocity is only extruded if collapsed
    135                 InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,true);
    136         }
    137 
    138130        //more output might be needed, when running in control
    139131        if(pKff0){
Note: See TracChangeset for help on using the changeset viewer.