|
| 1 | +#include "./MmemasstransportAnalysis.h" |
| 2 | +#include <math.h> |
| 3 | +#include "../toolkits/toolkits.h" |
| 4 | +#include "../classes/classes.h" |
| 5 | +#include "../classes/Inputs/TransientInput.h" |
| 6 | +#include "../classes/Inputs/TriaInput.h" |
| 7 | +#include "../classes/gauss/Gauss.h" |
| 8 | +#include "../shared/shared.h" |
| 9 | +#include "../modules/modules.h" |
| 10 | + |
| 11 | +/*Model processing*/ |
| 12 | +void MmemasstransportAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ |
| 13 | + /*No constraints*/ |
| 14 | +}/*}}}*/ |
| 15 | +void MmemasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ |
| 16 | + /*No loads*/ |
| 17 | +}/*}}}*/ |
| 18 | +void MmemasstransportAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ |
| 19 | + ::CreateNodes(nodes,iomodel,MmemasstransportAnalysisEnum,P1Enum); |
| 20 | +}/*}}}*/ |
| 21 | +int MmemasstransportAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ |
| 22 | + return 3; |
| 23 | +}/*}}}*/ |
| 24 | +void MmemasstransportAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ |
| 25 | + |
| 26 | + int nature=0; |
| 27 | + bool isdakota=0; |
| 28 | + |
| 29 | + /*Update elements: */ |
| 30 | + int counter=0; |
| 31 | + for(int i=0;i<iomodel->numberofelements;i++){ |
| 32 | + if(iomodel->my_elements[i]){ |
| 33 | + Element* element=(Element*)elements->GetObjectByOffset(counter); |
| 34 | + element->Update(inputs,i,iomodel,analysis_counter,analysis_type,P1Enum); |
| 35 | + counter++; |
| 36 | + } |
| 37 | + } |
| 38 | + |
| 39 | + /*Plug inputs into element:*/ |
| 40 | + iomodel->FetchDataToInput(inputs,elements,"md.geometry.thickness",ThicknessEnum); |
| 41 | + iomodel->FetchDataToInput(inputs,elements,"md.geometry.surface",SurfaceEnum); |
| 42 | + iomodel->FetchDataToInput(inputs,elements,"md.geometry.base",BaseEnum); |
| 43 | + iomodel->FetchDataToInput(inputs,elements,"md.initialization.sealevel",SealevelEnum,0); |
| 44 | + iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MmemasstransportMaskIceLevelsetEnum); |
| 45 | + iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MmemasstransportMaskOceanLevelsetEnum); |
| 46 | + iomodel->FetchDataToInput(inputs,elements,"md.mmemasstransport.thickness", MmemasstransportThicknessEnum); |
| 47 | + |
| 48 | + /*Initialize sea level cumulated sea level loads :*/ |
| 49 | + iomodel->ConstantToInput(inputs,elements,0.,AccumulatedDeltaIceThicknessEnum,P0Enum); |
| 50 | + iomodel->ConstantToInput(inputs,elements,0.,OldAccumulatedDeltaIceThicknessEnum,P0Enum); |
| 51 | + iomodel->ConstantToInput(inputs,elements,0.,DeltaIceThicknessEnum,P0Enum); |
| 52 | + |
| 53 | +}/*}}}*/ |
| 54 | +void MmemasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ |
| 55 | + |
| 56 | + int numoutputs; |
| 57 | + char** requestedoutputs = NULL; |
| 58 | + |
| 59 | + int nids,npart,nel; |
| 60 | + IssmDouble* ids=NULL; |
| 61 | + IssmDouble* partition = NULL; |
| 62 | + |
| 63 | + iomodel->FetchData(&nel,"md.mesh.numberofelements"); |
| 64 | + iomodel->FetchData(&ids,&nids,NULL,"md.mmemasstransport.ids"); |
| 65 | + //_printf_("nids: " << nids << "\n"); for(int i=0;i<nids;i++)_printf_(ids[i] << "|"); _printf_("\n"); |
| 66 | + parameters->AddObject(new DoubleMatParam(MmemasstransportModelidsEnum,ids,nids,1)); |
| 67 | + iomodel->FetchData(&partition,&npart,NULL,"md.mmemasstransport.partition"); |
| 68 | + if (npart!=nel)_error_("MmemasstransportAnalysis:UpdateParameters: partition vector should be distributed over elements, not vertices!"); |
| 69 | + parameters->AddObject(new DoubleMatParam(MmemasstransportPartitionEnum,partition,nel,1)); |
| 70 | + |
| 71 | + xDelete<IssmDouble>(ids); |
| 72 | + xDelete<IssmDouble>(partition); |
| 73 | + |
| 74 | + /*Requested outputs*/ |
| 75 | + iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.mmemasstransport.requested_outputs"); |
| 76 | + parameters->AddObject(new IntParam(MmemasstransportNumRequestedOutputsEnum,numoutputs)); |
| 77 | + if(numoutputs)parameters->AddObject(new StringArrayParam(MmemasstransportRequestedOutputsEnum,requestedoutputs,numoutputs)); |
| 78 | + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.mmemasstransport.requested_outputs"); |
| 79 | + |
| 80 | +}/*}}}*/ |
| 81 | + |
| 82 | +/*Finite Element Analysis*/ |
| 83 | +void MmemasstransportAnalysis::Core(FemModel* femmodel){/*{{{*/ |
| 84 | + _error_("not implemented"); |
| 85 | +}/*}}}*/ |
| 86 | +void MmemasstransportAnalysis::PreCore(FemModel* femmodel){/*{{{*/ |
| 87 | + _error_("not implemented"); |
| 88 | +}/*}}}*/ |
| 89 | +ElementVector* MmemasstransportAnalysis::CreateDVector(Element* element){/*{{{*/ |
| 90 | + /*Default, return NULL*/ |
| 91 | + return NULL; |
| 92 | +}/*}}}*/ |
| 93 | +ElementMatrix* MmemasstransportAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ |
| 94 | +_error_("Not implemented"); |
| 95 | +}/*}}}*/ |
| 96 | +ElementMatrix* MmemasstransportAnalysis::CreateKMatrix(Element* element){/*{{{*/ |
| 97 | + _error_("not implemented yet"); |
| 98 | +}/*}}}*/ |
| 99 | +ElementVector* MmemasstransportAnalysis::CreatePVector(Element* element){/*{{{*/ |
| 100 | +_error_("not implemented yet"); |
| 101 | +}/*}}}*/ |
| 102 | +void MmemasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/ |
| 103 | + |
| 104 | + /*do nothing:*/ |
| 105 | + return; |
| 106 | +}/*}}}*/ |
| 107 | +void MmemasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/ |
| 108 | + _error_("Not implemented yet"); |
| 109 | +}/*}}}*/ |
| 110 | +void MmemasstransportAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ |
| 111 | + |
| 112 | + /*Fetch number of nodes and dof for this finite element*/ |
| 113 | + IssmDouble time; |
| 114 | + IssmDouble ice[3]; |
| 115 | + IssmDouble ocean[3]; |
| 116 | + IssmDouble height; |
| 117 | + int numnodes = element->GetNumberOfNodes(); |
| 118 | + |
| 119 | + element->parameters->FindParam(&time,TimeEnum); |
| 120 | + |
| 121 | + TriaInput* h_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportThicknessEnum,time)); _assert_(h_input); |
| 122 | + TriaInput* i_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskIceLevelsetEnum,time)); _assert_(i_input); |
| 123 | + TriaInput* o_input=xDynamicCast<TriaInput*>(element->GetInput(MmemasstransportMaskOceanLevelsetEnum,time)); _assert_(o_input); |
| 124 | + |
| 125 | + Gauss* gauss=element->NewGauss(); |
| 126 | + for(int iv=0;iv<3;iv++){ |
| 127 | + gauss->GaussVertex(iv); |
| 128 | + i_input->GetInputValue(&ice[iv],gauss); |
| 129 | + o_input->GetInputValue(&ocean[iv],gauss); |
| 130 | + } |
| 131 | + h_input->GetInputAverage(&height); |
| 132 | + delete gauss; |
| 133 | + |
| 134 | + /*Add thickness ice and ocean levelsets as inputs to the tria element: */ |
| 135 | + element->AddInput(ThicknessEnum,&height,P0Enum); //very important, do not change the type of ThicknessEnum to P1 when it should be P0! |
| 136 | + element->AddInput(MaskIceLevelsetEnum,ice,P1Enum); |
| 137 | + element->AddInput(MaskOceanLevelsetEnum,ocean,P1Enum); |
| 138 | + |
| 139 | +}/*}}}*/ |
| 140 | +void MmemasstransportAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ |
| 141 | + /*Default, do nothing*/ |
| 142 | + return; |
| 143 | +}/*}}}*/ |
0 commit comments