source: issm/trunk-jpl/src/modules/InternalFront/InternalFront.cpp@ 13236

/*\file InternalFront.c
 *\brief: build pressureload
 */

#include "./InternalFront.h"

void InternalFrontUsage(void) {/*{{{*/
        _pprintLine_("");
        _pprintLine_("   usage: icefront = " << __FUNCT__ << "(md);");
        _pprintLine_("");
}/*}}}*/
WRAPPER(InternalFront){

        bool*   elementonwater=NULL;
        int*    elements=NULL;
        int*    connectivity=NULL;
        int*    elementconnectivity=NULL;
        int*    front=NULL;
        double* front2=NULL;
        bool    found;
        int     numberofelements,numberofsegments;
        int     N,M;
        int     i,j,ii,jj,id;

        /*Boot module: */
        MODULEBOOT();

        /*checks on arguments on the matlab side: */
        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InternalFrontUsage);

        /*Fetch required fields*/
        FetchData(&numberofelements,mxGetAssignedField(MODEL,0,"numberofelements"));
        if(numberofelements<=0) _error_("No elements found in the model");
        FetchData(&elements,&M,&N,mxGetAssignedField(MODEL,0,"elements"));
        if(M!=numberofelements || N!=3) _error_("Field 'elements' should be of size [md.numberofelements 3]");
        FetchData(&elementonwater,&M,&N,mxGetAssignedField(MODEL,0,"elementonwater"));
        if(M!=numberofelements || N!=1) _error_("Field 'elementonwater' should be of size [md.numberofelements 1]");
        FetchData(&elementconnectivity,&M,&N,mxGetAssignedField(MODEL,0,"elementconnectivity"));
        if(M!=numberofelements || N!=3) _error_("Field 'elementconnectivity' should be of size [md.numberofelements 3]");

        /*Allocate and initialize all variables*/
        numberofsegments=0;
        front=xNew<int>(3*numberofelements*4);

        /*Loop over all elements on water*/
        for(i=0;i<numberofelements;i++){

                /*Skip if on water*/
                if(!elementonwater[i]) continue;

                /*Loop over all three adjacent elements*/
                for(j=0;j<3;j++){

                        /*Skip if adjacent element does not exist or is on water*/
                        id=elementconnectivity[i*3+j];
                        if(id==0) continue;
                        if(elementonwater[id-1]) continue;

                        /*We have an ice front to add here, let's go!*/
                        for(ii=0;ii<3;ii++){

                                found=false;
                                for(jj=0;jj<3;jj++){
                                        if(elements[(id-1)*3+ii]==elements[i*3+jj]){
                                                found=true;
                                                break;
                                        }
                                }

                                /*OK, we just found the node of id, which is not in i. We have the segment*/
                                if(!found){
                                        front[numberofsegments*4+0]=elements[(id-1)*3+(ii+1)%3];
                                        front[numberofsegments*4+1]=elements[(id-1)*3+(ii+2)%3];
                                        front[numberofsegments*4+2]=id;
                                        front[numberofsegments*4+3]=IceEnum;
                                        numberofsegments+=1;
                                        break;
                                }
                        }

                        /*In debugging mode, check that the segment has been found*/
                        _assert_(!found);
                }
        }

        /*Now that we know how many segments there is we can allocate the final matrix*/
        if(numberofsegments){
                front2=xNew<double>(4*numberofsegments);
                for(i=0;i<4*numberofsegments;i++) front2[i]=(double)front[i];
        }
        xfree((void**)&front);

        /*write output datasets: */
        WriteData(FRONT,front2,numberofsegments,4);

        /*end module: */
        MODULEEND();
}