Changeset 20665


Ignore:
Timestamp:
05/29/16 17:46:38 (9 years ago)
Author:
felicity
Message:

CHG: using diffusion equation to extrude rather than advection

Location:
issm/trunk-jpl/src/c
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp

    r18929 r20665  
    4949ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
    5050
    51         /*compute all stiffness matrices for this element*/
    52         ElementMatrix* Ke1=CreateKMatrixVolume(element);
    53         ElementMatrix* Ke2=CreateKMatrixSurface(element);
    54         ElementMatrix* Ke3=CreateKMatrixBed(element);
    55         ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    56 
    57         /*clean-up and return*/
    58         delete Ke1;
    59         delete Ke2;
    60         delete Ke3;
    61         return Ke;
    62 }/*}}}*/
    63 ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
    64 
    65         if(!element->IsOnBase()) return NULL;
    66 
    67         /*Intermediaries */
    68         int         dim;
    69         IssmDouble  Jdet,D,normal[3];
    70         IssmDouble *xyz_list_base = NULL;
    71 
    72         /*Get dimension*/
    73         element->FindParam(&dim,DomainDimensionEnum);
    74 
    75         /*Fetch number of nodes and dof for this finite element*/
    76         int numnodes = element->GetNumberOfNodes();
    77 
    78         /*Initialize Element vector and other vectors*/
    79         ElementMatrix* Ke    = element->NewElementMatrix();
    80         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    81 
    82         /*Retrieve all inputs and parameters*/
    83         element->GetVerticesCoordinatesBase(&xyz_list_base);
    84         element->NormalBase(&normal[0],xyz_list_base);
    85 
    86         /* Start  looping on the number of gaussian points: */
    87         Gauss* gauss=element->NewGaussBase(2);
    88         for(int ig=gauss->begin();ig<gauss->end();ig++){
    89                 gauss->GaussPoint(ig);
    90 
    91                 element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
    92                 element->NodalFunctions(basis,gauss);
    93                 D = - gauss->weight*Jdet*normal[dim-1];
    94 
    95                 TripleMultiply(basis,1,numnodes,1,
    96                                         &D,1,1,0,
    97                                         basis,1,numnodes,0,
    98                                         &Ke->values[0],1);
    99         }
    100 
    101         /*Clean up and return*/
    102         xDelete<IssmDouble>(xyz_list_base);
    103         xDelete<IssmDouble>(basis);
    104         delete gauss;
    105         return Ke;
    106 }
    107 /*}}}*/
    108 ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
    109 
    110         if(!element->IsOnSurface()) return NULL;
    111 
    112         /*Intermediaries */
    113         int         dim;
    114         IssmDouble  Jdet,D,normal[3];
    115         IssmDouble *xyz_list_top = NULL;
    116 
    117         /*Get dimension*/
    118         element->FindParam(&dim,DomainDimensionEnum);
    119 
    120         /*Fetch number of nodes and dof for this finite element*/
    121         int numnodes = element->GetNumberOfNodes();
    122 
    123         /*Initialize Element vector and other vectors*/
    124         ElementMatrix* Ke    = element->NewElementMatrix();
    125         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    126 
    127         /*Retrieve all inputs and parameters*/
    128         element->GetVerticesCoordinatesTop(&xyz_list_top);
    129         element->NormalTop(&normal[0],xyz_list_top);
    130 
    131         /* Start  looping on the number of gaussian points: */
    132         Gauss* gauss=element->NewGaussTop(2);
    133         for(int ig=gauss->begin();ig<gauss->end();ig++){
    134                 gauss->GaussPoint(ig);
    135 
    136                 element->JacobianDeterminantTop(&Jdet,xyz_list_top,gauss);
    137                 element->NodalFunctions(basis,gauss);
    138                 D = - gauss->weight*Jdet*normal[dim-1];
    139 
    140                 TripleMultiply(basis,1,numnodes,1,
    141                                         &D,1,1,0,
    142                                         basis,1,numnodes,0,
    143                                         &Ke->values[0],1);
    144         }
    145 
    146         /*Clean up and return*/
    147         xDelete<IssmDouble>(xyz_list_top);
    148         xDelete<IssmDouble>(basis);
    149         delete gauss;
    150         return Ke;
    151 }
    152 /*}}}*/
    153 ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
    154 
    15551        /*Intermediaries */
    15652        IssmDouble  Jdet,D;
    … …  
    16662        /*Initialize Element vector and other vectors*/
    16763        ElementMatrix* Ke     = element->NewElementMatrix();
    168         IssmDouble*    B      = xNew<IssmDouble>(numnodes);
    169         IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
     64        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    17065
    17166        /*Retrieve all inputs and parameters*/
    … …  
    17873
    17974                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    180                 element->NodalFunctions(Bprime,gauss);
    181                 GetB(B,element,dim,xyz_list,gauss);
    18275                D=gauss->weight*Jdet;
     76                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    18377
    184                 TripleMultiply(B,1,numnodes,1,
    185                                         &D,1,1,0,
    186                                         Bprime,1,numnodes,0,
    187                                         &Ke->values[0],1);
     78                for(int i=0;i<numnodes;i++){
     79                        for(int j=0;j<numnodes;j++){
     80                                Ke->values[i*numnodes+j] += gauss->weight*Jdet*(dbasis[(dim-1)*numnodes+i]*dbasis[(dim-1)*numnodes+j]);
     81                        }
     82                }
    18883        }
    18984
    19085        /*Clean up and return*/
    19186        xDelete<IssmDouble>(xyz_list);
    192         xDelete<IssmDouble>(B);
    193         xDelete<IssmDouble>(Bprime);
    19487        delete gauss;
    19588        return Ke;
    196 }
    197 /*}}}*/
     89}/*}}}*/
    19890ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/
    19991        return NULL;
    20092}/*}}}*/
    201 void           ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    202         /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
    203                 where hi is the interpolation function for node i.*/
    204 
    205         /*Fetch number of nodes for this finite element*/
    206         int numnodes = element->GetNumberOfNodes();
    207 
    208         /*Get nodal functions derivatives*/
    209         IssmDouble* dbasis=xNew<IssmDouble>(dim*numnodes);
    210         element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    211 
    212         /*Build B: */
    213         for(int i=0;i<numnodes;i++){
    214                 B[i] = dbasis[(dim-1)*numnodes+i];
    215         }
    216 
    217         /*Clean-up*/
    218         xDelete<IssmDouble>(dbasis);
    219 }
    220 /*}}}*/
    22193void           ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
    22294           _error_("not implemented yet");

  • issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h

    r18929 r20665  
    2525                ElementMatrix* CreateJacobianMatrix(Element* element);
    2626                ElementMatrix* CreateKMatrix(Element* element);
    27                 ElementMatrix* CreateKMatrixBed(Element* element);
    28                 ElementMatrix* CreateKMatrixSurface(Element* element);
    29                 ElementMatrix* CreateKMatrixVolume(Element* element);
    3027                ElementVector* CreatePVector(Element* element);
    31                 void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    3228                void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
    3329                void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r20645 r20665  
    642642        /*Get viscosity*/
    643643        this->GetViscosity(&viscosity,eps_eff);
     644        _assert_(!xIsNan<IssmDouble>(viscosity));
    644645
    645646        /*Assign output pointer*/

  • issm/trunk-jpl/src/c/classes/Vertex.cpp

    r19984 r20665  
    164164                        this->y = newy;
    165165                        vy->SetValue(this->pid,vely,INS_VAL);
     166                        _assert_(!xIsNan<IssmDouble>(vely));
    166167                        return;
    167168                case Domain3DEnum:
    … …  
    171172                        this->z = newz;
    172173                        vz->SetValue(this->pid,velz,INS_VAL);
     174                        _assert_(!xIsNan<IssmDouble>(velz));
    173175                        return;
    174176                default:
Note: See TracChangeset for help on using the changeset viewer.