source: issm/branches/trunk-larour-SLPS2022/src/c/classes/GrdLoads.cpp@ 27104

/*
 * \file GrdLoads.cpp
 * \brief: Implementation of GrdLoads class
 */

/*Headers: {{{*/
#ifdef HAVE_CONFIG_H
        #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./GrdLoads.h"
#include "./SealevelGeometry.h"
using namespace std;
/*}}}*/

/*Object constructors and destructor*/
GrdLoads::GrdLoads(int nel,SealevelGeometry* slgeom){ /*{{{*/

        /*allocate:*/
        vloads=new Vector<IssmDouble>(nel);
        for (int i=0;i<SLGEOM_NUMLOADS;i++) vsubloads[i]=new Vector<IssmDouble>(slgeom->nbar[i]);

        vsealevelloads=new Vector<IssmDouble>(nel);
        vsealevelloads->Set(0);vsealevelloads->Assemble();

        vsubsealevelloads=new Vector<IssmDouble>(slgeom->nbar[SLGEOM_OCEAN]);

        /*make sure all pointers that are not allocated are NULL:*/
        loads=NULL;
        for (int i=0;i<SLGEOM_NUMLOADS;i++)subloads[i]=NULL;
        sealevelloads=NULL;
        subsealevelloads=NULL;


}; /*}}}*/
GrdLoads::~GrdLoads(){ /*{{{*/

        delete vloads;
        xDelete<IssmDouble>(loads);
        for(int i=0;i<SLGEOM_NUMLOADS;i++){
                delete vsubloads[i];
                xDelete<IssmDouble>(subloads[i]);
        }
        delete vsealevelloads;
        xDelete<IssmDouble>(sealevelloads);
        delete vsubsealevelloads;
        xDelete<IssmDouble>(subsealevelloads);
}; /*}}}*/

void GrdLoads::BroadcastLoads(void){ /*{{{*/

        /*Initialize barycentre vectors, now that we know their size: */
        vloads->Assemble();
        for (int i=0;i<SLGEOM_NUMLOADS;i++){
                vsubloads[i]->Assemble();
        }

        /*Avoid leaks:*/
        if(loads)xDelete<IssmDouble>(loads);
        for (int i=0;i<SLGEOM_NUMLOADS;i++){
                if(subloads[i])xDelete<IssmDouble>(subloads[i]);
        }

        /*Serialize:*/
        loads=vloads->ToMPISerial();
        for (int i=0;i<SLGEOM_NUMLOADS;i++){
                subloads[i]=vsubloads[i]->ToMPISerial();
        }

} /*}}}*/
void GrdLoads::AssembleSealevelLoads(void){ /*{{{*/

        vsealevelloads->Assemble();
        vsubsealevelloads->Assemble();

} /*}}}*/
void GrdLoads::BroadcastSealevelLoads(void){ /*{{{*/

        /*Avoid leakds:*/
        if(sealevelloads)xDelete<IssmDouble>(sealevelloads);
        if(subsealevelloads)xDelete<IssmDouble>(subsealevelloads);

        /*Serialize:*/
        sealevelloads=vsealevelloads->ToMPISerial();
        subsealevelloads=vsubsealevelloads->ToMPISerial();

} /*}}}*/
void GrdLoads::SHDegree2Coefficients(IssmDouble* deg2coeff, FemModel* femmodel, SealevelGeometry* slgeom){ /*{{{*/

        IssmDouble area,re, S;
        int ylmindex, intj;
        IssmDouble deg2coeff_local[5];

        femmodel->parameters->FindParam(&re,SolidearthPlanetRadiusEnum);

        for (int c=0;c<5;c++) deg2coeff_local[c]=0;

        for(Object* & object : femmodel->elements->objects){
                Element* element = xDynamicCast<Element*>(object);
                //first, loads on centroids
                S=0;
                S+=loads[element->Sid()];

                if(sealevelloads) S+=sealevelloads[element->Sid()];
                if(S!=0){
                        element->Element::GetInputValue(&area,AreaEnum);

                        for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients: 2,0; 2,1cos; 2,1sin; 2,2cos; 2,2sin
                                ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
                                deg2coeff_local[c] += S*area/re/re*slgeom->Ylm[ylmindex];
                        }
                }
                //add loads on subelement barycenters
                for (int i=0;i<SLGEOM_NUMLOADS;i++){
                        if (slgeom->issubelement[i][element->lid]){
                                intj=slgeom->subelementmapping[i][element->lid];
                                S=0;
                                S+=subloads[i][intj];
                                if(i==SLGEOM_OCEAN && sealevelloads) S+=subsealevelloads[intj];
                                if(S!=0){
                                        area=slgeom->area_subel[i][intj];
                                        for (int c=0;c<5;c++){ //degree l=2 has 2*l+1=5 coefficients
                                                ylmindex=(4+c)*slgeom->localnel+element->lid; // starting at index=l^2
                                                deg2coeff_local[c] += S*area/re/re*slgeom->Ylm_subel[i][ylmindex];
                                        }
                                }
                        }
                }
        }

        //sum each degree 2 coefficient across all cpus to get global total
        ISSM_MPI_Reduce (&deg2coeff_local[0],&deg2coeff[0],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
        ISSM_MPI_Bcast(&deg2coeff[0],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
        ISSM_MPI_Reduce (&deg2coeff_local[1],&deg2coeff[1],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
        ISSM_MPI_Bcast(&deg2coeff[1],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
        ISSM_MPI_Reduce (&deg2coeff_local[2],&deg2coeff[2],1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() );
        ISSM_MPI_Bcast(&deg2coeff[2],1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());

} /*}}}*/