Ignore:
Timestamp:
04/20/10 16:35:48 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added discontinuous Galerkin for Balanced thickness solution

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp

    r3567 r3582  
    88#include "../../shared/shared.h"
    99#include "../../include/macros.h"
     10#include "../../include/typedefs.h"
    1011#include "../IoModel.h"
    11 
    1212
    1313void    CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){
    1414
    15         DataSet*    loads    = NULL;
     15        /*Intermediary*/
     16        int i,j;
     17        int i1,i2;
     18        int pos1,pos2;
     19        double e1,e2;
     20
     21        /*Output*/
     22        DataSet* loads=NULL;
     23
     24        /*numericalflux intermediary data: */
     25        char numericalflux_type[NUMERICALFLUXSTRING];
     26        int  numericalflux_id;
     27        int  numericalflux_node_ids[MAX_NUMERICALFLUX_NODES];
     28        int  numericalflux_elem_id;
     29        double numericalflux_h[MAX_NUMERICALFLUX_NODES];
    1630
    1731        /*Create loads: */
    1832        loads   = new DataSet(LoadsEnum);
    19        
    20        
     33
     34        /*Get edges and elements*/
     35        IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges");
     36        IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");
     37        IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");
     38
     39        /*First load data:*/
     40        for (i=0;i<iomodel->numberofedges;i++){
     41
     42                /*Get left and right elements*/
     43                e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2]
     44                e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2]
     45
     46                /*Now, if this element is not in the partition, pass: */
     47                if(!iomodel->my_elements[(int)e1]) continue;
     48
     49                /*Create load*/
     50                numericalflux_id=i+1; //Matlab indexing
     51                numericalflux_elem_id=(int)e1+1;//id is in matlab index
     52
     53                /*1: Get vertices ids*/
     54                i1=(int)iomodel->edges[4*i+0];
     55                i2=(int)iomodel->edges[4*i+1];
     56
     57                if (!isnan(e2)){
     58                        strcpy(numericalflux_type,"internal");
     59
     60                        /*Now, we must get the nodes of the 4 nodes located on the edge*/
     61
     62                        /*2: Get the column where these ids are located in the index*/
     63                        pos1=pos2=UNDEF;
     64                        for(j=0;j<3;j++){
     65                                if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
     66                                if (iomodel->elements[3*(int)e2+j]==i1) pos2=j+1;
     67                        }
     68                        ISSMASSERT(pos1!=UNDEF && pos2!=UNDEF);
     69
     70                        /*3: We have the id of the elements and the position of the vertices in the index
     71                         * we can compute their dofs!*/
     72                        numericalflux_node_ids[0]=3*(int)e1+pos1;       //ex: 1 2 3
     73                        numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1; //ex: 2 3 1
     74                        numericalflux_node_ids[2]=3*(int)e2+pos2;           //ex: 1 2 3
     75                        numericalflux_node_ids[3]=3*(int)e2+((pos2+1)%3)+1; //ex: 3 1 2
     76
     77                        numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1] -1];
     78                        numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
     79                        numericalflux_h[2]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[2]-1]-1];
     80                        numericalflux_h[3]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[3]-1]-1];
     81                }
     82                else{
     83                        strcpy(numericalflux_type,"boundary");
     84
     85                        /*2: Get the column where these ids are located in the index*/
     86                        pos1==UNDEF;
     87                        for(j=0;j<3;j++){
     88                                if (iomodel->elements[3*(int)e1+j]==i1) pos1=j+1;
     89                        }
     90                        ISSMASSERT(pos1!=UNDEF);
     91
     92                        /*3: We have the id of the elements and the position of the vertices in the index
     93                         * we can compute their dofs!*/
     94                        numericalflux_node_ids[0]=3*(int)e1+pos1;
     95                        numericalflux_node_ids[1]=3*(int)e1+(pos1%3)+1;
     96
     97                        numericalflux_h[0]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[0]-1]-1];
     98                        numericalflux_h[1]=iomodel->thickness[(int)iomodel->elements[numericalflux_node_ids[1]-1]-1];
     99                        numericalflux_h[2]=UNDEF;
     100                        numericalflux_h[3]=UNDEF;
     101                }
     102
     103                loads->AddObject(new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id,numericalflux_h));
     104        }
     105
     106        /*Free data: */
     107        xfree((void**)&iomodel->edges);
     108        xfree((void**)&iomodel->elements);
     109        xfree((void**)&iomodel->thickness);
     110
    21111        /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these
    22112         * datasets, it will not be redone: */
    23113        loads->Presort();
    24 
    25         cleanup_and_return:
    26114
    27115        /*Assign output pointer: */
    … …  
    29117
    30118}
    31 
    32 
Note: See TracChangeset for help on using the changeset viewer.