source: issm/trunk-jpl/src/c/modules/ModelProcessorx/StressbalanceSIA/CreateConstraintsStressbalanceSIA.cpp@ 16521

/*
 * CreateConstraintsStressbalanceSIA.c:
 */

#include "../../../toolkits/toolkits.h"
#include "../../../classes/classes.h"
#include "../../../shared/shared.h"
#include "../ModelProcessorx.h"

void    CreateConstraintsStressbalanceSIA(Constraints** pconstraints, IoModel* iomodel){

        /*Intermediary*/
        int        count;
        IssmDouble yts;
        bool       isSIA;

        /*Output*/
        Constraints* constraints = NULL;

        /*Recover pointer: */
        constraints=*pconstraints;

        /*Fetch parameters: */
        iomodel->Constant(&yts,ConstantsYtsEnum);
        iomodel->Constant(&isSIA,FlowequationIsSIAEnum);

        /*Now, is the flag isSIA on? otherwise, do nothing: */
        if (!isSIA) return;

        /*Fetch data: */
        iomodel->FetchData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);

        /*Initialize conunter*/
        count=0;

        /*vx and vy are spc'd if we are not on nodeonSIA: */
        for(int i=0;i<iomodel->numberofvertices;i++){
                /*keep only this partition's nodes:*/
                if((iomodel->my_vertices[i])){
                        if (reCast<int,IssmDouble>(iomodel->Data(FlowequationVertexEquationEnum)[i])!=SIAApproximationEnum){

                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,0,StressbalanceSIAAnalysisEnum));
                                count++;

                                constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,0,StressbalanceSIAAnalysisEnum));
                                count++;
                        }
                        else{
                                if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvxEnum)[i])){
                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->Data(StressbalanceSpcvxEnum)[i]/yts,StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 1 to vx.
                                        count++;
                                }

                                if (!xIsNan<IssmDouble>(iomodel->Data(StressbalanceSpcvyEnum)[i])){
                                        constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,2,iomodel->Data(StressbalanceSpcvyEnum)[i]/yts,StressbalanceSIAAnalysisEnum)); //add count'th spc, on node i+1, setting dof 2 to vy
                                        count++;
                                }
                        }
                }
        }

        /*Free data: */
        iomodel->DeleteData(3,StressbalanceSpcvxEnum,StressbalanceSpcvyEnum,FlowequationVertexEquationEnum);

        /*Assign output pointer: */
        *pconstraints=constraints;
}