Changeset 5720

Timestamp:

09/09/10 12:58:29 (15 years ago)

Author:

Mathieu Morlighem

Message:

Spawn tria for MacAyeal 3d ice front

File:

• issm/trunk/src/c/objects/Loads/Icefront.cpp (modified) (1 diff)

issm/trunk/src/c/objects/Loads/Icefront.cpp

@@ -504 +504 @@
 void Icefront::CreatePVectorDiagnosticMacAyeal3d( Vec pg){
 
-        int       i;
-        const int numnodes= NUMVERTICESSEG;
-        const int numdofs = numnodes *NDOF2;
-        int*      doflist=NULL;
-        double    xyz_list[numnodes][3];
-
-        /*Objects: */
-        double    pe_g[numdofs] = {0.0};
-
-        /*Segment*/
-        double   normal[2];
-        int      num_gauss;
-        double*  segment_gauss_coord=NULL;
-        double*  gauss_weights=NULL;
-        double   gauss_weight;
-        double   gauss_coord;
-        int      ig;
-        double   Jdet;
-        double   L[2];
-        double   thickness,bed,pressure;
-        double   ice_pressure,water_pressure,air_pressure,rho_water,rho_ice,gravity;
-        double   surface_under_water,base_under_water;
-        int      fill;
-
-        /*Inputs*/
-        Input*  thickness_input=NULL;
-        Input*  bed_input=NULL;
-        Penta*  penta=NULL;
-
-        BeamRef* beam=NULL;
-
-        //check that the element is onbed (collapsed formulation) otherwise:pe=0
-        if(!element->GetOnBed()) return;
-
-        /*Recover inputs: */
-        thickness_input=((Penta*)element)->inputs->GetInput(ThicknessEnum);
-        bed_input      =((Penta*)element)->inputs->GetInput(BedEnum);
-        inputs->GetParameterValue(&fill,FillEnum);
-
-        /* Get dof list and node coordinates: */
-        GetDofList(&doflist,MacAyealApproximationEnum);
-        GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes);
-
-        /*Get Matpar parameters*/
-        rho_water=matpar->GetRhoWater();
-        rho_ice  =matpar->GetRhoIce();
-        gravity  =matpar->GetG();
-
-        /*Get normal*/
-        GetNormal(&normal[0],xyz_list);
-
-        //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions
-        num_gauss=3;
-        gaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss);
-
-        for(ig=0;ig<num_gauss;ig++){
-
-                /*Pick up the gaussian point: */
-                gauss_weight=gauss_weights[ig];
-                gauss_coord=segment_gauss_coord[ig];
-
-                //compute Jacobian of segment
-                beam->GetJacobianDeterminant2d(&Jdet,&xyz_list[0][0],gauss_coord);
-
-                /*Get thickness and bed*/
-                this->GetParameterValue(&thickness,gauss_coord,thickness_input);
-                this->GetParameterValue(&bed,      gauss_coord,bed_input);
-
-                /*Compute total pressure applied to ice front*/
-                if (fill==WaterEnum){
-                        //icefront ends in water: 
-                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-                        air_pressure=0;
-
-                        //Now deal with water pressure
-                        surface_under_water=min(0,thickness+bed); // 0 if the top of the glacier is above water level
-                        base_under_water=min(0,bed);              // 0 if the bottom of the glacier is above water level
-                        water_pressure=1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
-                }
-                else if (fill==AirEnum){
-                        ice_pressure=1.0/2.0*gravity*rho_ice*pow(thickness,2);
-                        air_pressure=0;
-                        water_pressure=0;
-                }
-                else{
-                        ISSMERROR("fill type %i not supported yet",fill);
-                }
-                pressure = ice_pressure + water_pressure + air_pressure;
-
-                /*Get L*/
-                beam->GetNodalFunctions(&L[0],gauss_coord);
-
-                for (i=0;i<numnodes;i++){
-                        pe_g[2*i+0]+= pressure*Jdet*gauss_weights[ig]*normal[0]*L[i];
-                        pe_g[2*i+1]+= pressure*Jdet*gauss_weights[ig]*normal[1]*L[i];
-                }
-        } //for(ig=0;ig<num_gauss;ig++)
-
-
-        /*Plug pe_g into vector: */
-        VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES);
-
-        /*Free ressources:*/
-        xfree((void**)&segment_gauss_coord);
-        xfree((void**)&gauss_weights);
-        xfree((void**)&doflist);
-
+        Icefront* icefront=NULL;
+        Penta*    penta=NULL;
+        Tria*     tria=NULL;
+        bool      onbed;
+
+        /*Cast element onto Penta*/
+        penta   =(Penta*)this->element;
+
+        /*Return if not on bed*/
+        penta->inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+        if(!onbed) return;
+
+        /*Spawn Tria and call MacAyeal2d*/
+        tria    =(Tria*)penta->SpawnTria(0,1,2);
+        icefront=(Icefront*)this->copy();
+        icefront->element=tria;
+        icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum));
+        icefront->CreatePVectorDiagnosticMacAyeal2d(pg);
+
+        /*clean-up*/
+        delete tria->matice;
+        delete tria;
+        delete icefront;
 }
 /*}}}*/