Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24710) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 24711) @@ -1006,4 +1006,138 @@ /*Update Levelset*/ this->AddInput2(distanceenum,&ls[0],P1Enum); +} +/*}}}*/ +void Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/ + + _assert_(end_time>start_time); + + /*Intermediaries*/ + IssmDouble averaged_values[NUMVERTICES]; + IssmDouble current_values[NUMVERTICES]; + IssmDouble dt; + int found,start_offset,end_offset; + int averaging_method=0; + + + /*Get transient input time steps*/ + int numtimesteps; + IssmDouble *timesteps = NULL; + TransientInput2* transient_input = this->inputs2->GetTransientInput(transientinput_enum); + transient_input->GetAllTimes(×teps,&numtimesteps); + + /*go through the timesteps, and grab offset for start and end*/ + found=binary_search(&start_offset,start_time,timesteps,numtimesteps); + if(!found) _error_("Input not found (is TransientInput sorted ?)"); + found=binary_search(&end_offset,end_time,timesteps,numtimesteps); + if(!found) _error_("Input not found (is TransientInput sorted ?)"); + + Gauss* gauss=this->NewGauss(); + + /*stack the input for each timestep in the slice*/ + int offset = start_offset; + while(offset < end_offset ){ + + + if(offset==-1){ + /*get values for the first time: */ + _assert_(start_timeGetPentaInput(0); + + int interpolation = input->GetInterpolation(); + _assert_(interpolation==P1Enum); + /*Intermediaries*/ + int numindices; + int indices[NUMVERTICES]; + numindices = NUMVERTICES; + for(int i=0;ilid; + input->Serve(numindices,&indices[0]); + + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(¤t_values[iv],gauss); + } + dt = timesteps[0] - start_time; _assert_(dt>0.); + } + else{ + PentaInput2* input = transient_input->GetPentaInput(offset+1); + int interpolation = input->GetInterpolation(); + _assert_(interpolation==P1Enum); + /*Intermediaries*/ + int numindices; + int indices[NUMVERTICES]; + numindices = NUMVERTICES; + for(int i=0;ilid; + input->Serve(numindices,&indices[0]); + + for(int iv=0;ivGaussVertex(iv); + input->GetInputValue(¤t_values[iv],gauss); + } + if(offset == numtimesteps-1){ + dt = end_time - timesteps[offset]; _assert_(dt>0.); + } + else{ + dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.); + } + } + + switch(averaging_method){ + case 0: /*Arithmetic mean*/ + if(offset==start_offset){ + for(int iv=0;iv1.e-50); + averaged_values[iv] = dt*1./current_values[iv]; + } + } + else{ + for(int iv=0;iv1.e-50); + averaged_values[iv] += dt*1./current_values[iv]; + } + } + break; + default: + _error_("averaging method is not recognised"); + } + + offset+=1; + } + + /*Integration done, now normalize*/ + switch(averaging_method){ + case 0: //Arithmetic mean + for(int iv=0;ivAddInput2(averagedinput_enum,&averaged_values[0],P1Enum); + + /*Cleanup*/ + delete gauss; + xDelete(timesteps); } /*}}}*/ @@ -2173,5 +2307,5 @@ delete gauss; return flux; - + } /*}}}*/ @@ -2929,5 +3063,5 @@ /*First, serarch the input: */ - Input2* data=this->GetInput2(natureofdataenum); + Input2* data=this->GetInput2(natureofdataenum); /*figure out if we have the vertex id: */ @@ -4041,9 +4175,9 @@ /*Get material parameters :*/ rho_ice=FindParam(MaterialsRhoIceEnum); - Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input); + Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input); Input2* gllevelset_input = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input); Input2* scalefactor_input = NULL; if(scaled==true){ - scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); + scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); } ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -4090,5 +4224,5 @@ Input2* scalefactor_input = NULL; if(scaled==true){ - scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); + scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); } ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24710) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 24711) @@ -68,4 +68,5 @@ void CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum); ElementMatrix* CreateBasalMassMatrix(void); + void CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time); void ElementResponse(IssmDouble* presponse,int response_enum); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); Index: /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 24710) +++ /issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp (revision 24711) @@ -138,5 +138,4 @@ void TransientInput2::AddTriaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/ - /*Check whether this is the last time step that we have*/ if(this->numtimesteps){ @@ -181,5 +180,41 @@ void TransientInput2::AddPentaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/ - _error_("not implemented yet, look at TransientInput2::AddTriaTimeInput"); + /*Check whether this is the last time step that we have*/ + if(this->numtimesteps){ + if(fabs(this->timesteps[this->numtimesteps-1]-time)<1.0e-5){ + this->AddPentaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in); + return; + } + } + + /*This is a new time step! we need to add it to the list*/ + if(this->numtimesteps>0 && timetimesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially"); + + IssmDouble *old_timesteps = NULL; + Input2 **old_inputs = NULL; + if (this->numtimesteps > 0){ + old_timesteps=xNew(this->numtimesteps); + xMemCpy(old_timesteps,this->timesteps,this->numtimesteps); + xDelete(this->timesteps); + old_inputs=xNew(this->numtimesteps); + xMemCpy(old_inputs,this->inputs,this->numtimesteps); + xDelete(this->inputs); + } + + this->numtimesteps=this->numtimesteps+1; + this->timesteps=xNew(this->numtimesteps); + this->inputs = xNew(this->numtimesteps); + + if (this->numtimesteps > 1){ + xMemCpy(this->inputs,old_inputs,this->numtimesteps-1); + xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1); + xDelete(old_timesteps); + xDelete(old_inputs); + } + + /*go ahead and plug: */ + this->timesteps[this->numtimesteps-1] = time; + this->inputs[this->numtimesteps-1] = NULL; + this->AddPentaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in); } Index: /issm/trunk-jpl/test/NightlyRun/runme.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/runme.py (revision 24710) +++ /issm/trunk-jpl/test/NightlyRun/runme.py (revision 24711) @@ -157,5 +157,5 @@ print(("File {} saved. \n".format(os.path.join('..', 'Archives', archive_name + '.arch')))) - #ELSE: CHECK TEST + #ELSE: CHECK TEST else: #load archive Index: /issm/trunk-jpl/test/NightlyRun/test701.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test701.py (revision 24710) +++ /issm/trunk-jpl/test/NightlyRun/test701.py (revision 24711) @@ -13,6 +13,4 @@ md = bamgflowband(model(), x, b + h, b, 'hmax', 80.) -print(isinstance(md, model)) - #Geometry #interp1d returns a function to be called on md.mesh.x md.geometry.surface = np.interp(md.mesh.x, x, b + h) @@ -21,18 +19,18 @@ #mask -md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices, )) +md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices)) md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0. -md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices, )) - 0.5 +md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices)) - 0.5 #materials -md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, )) +md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices)) md.materials.rheology_B = paterson(md.initialization.temperature) -md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, )) +md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements)) #friction -md.friction.coefficient = np.zeros((md.mesh.numberofvertices, )) +md.friction.coefficient = np.zeros((md.mesh.numberofvertices)) md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20. -md.friction.p = np.ones((md.mesh.numberofelements, )) -md.friction.q = np.ones((md.mesh.numberofelements, )) +md.friction.p = np.ones((md.mesh.numberofelements)) +md.friction.q = np.ones((md.mesh.numberofelements)) #Boundary conditions @@ -47,6 +45,4 @@ #Misc -print(isinstance(md, model)) -print(type(md)) md = setflowequation(md, 'FS', 'all') md.stressbalance.abstol = np.nan @@ -67,5 +63,5 @@ for i in ['MINI', 'MINIcondensed', 'TaylorHood', 'LATaylorHood', 'CrouzeixRaviart', 'LACrouzeixRaviart']: print(' ') - print(' == == == Testing ' + i + ' Full - Stokes Finite element == == = ') + print('=====Testing ' + i + ' Full-Stokes Finite element=====') md.flowequation.fe_FS = i md = solve(md, 'Stressbalance')