Changeset 15744


Ignore:
Timestamp:
08/07/13 15:17:22 (12 years ago)
Author:
seroussi
Message:

CHG: trying to fix coupling HO/FS

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Penta.cpp

    r15741 r15744  
    401401
    402402        if(Ke){
    403                 /*Condense if requested*/
    404                 if(this->element_type==MINIcondensedEnum){
    405                         int indices[3]={18,19,20};
    406                         Ke->StaticCondensation(3,&indices[0]);
    407                 }
    408                 else if(this->element_type==P1bubblecondensedEnum){
    409                         int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    410                         int offset = 0;
    411                         for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    412                         int* indices=xNew<int>(size);
    413                         for(int i=0;i<size;i++) indices[i] = offset+i;
    414                         Ke->StaticCondensation(size,indices);
    415                         xDelete<int>(indices);
     403                int approximation;
     404                inputs->GetInputValue(&approximation,ApproximationEnum);
     405                if(approximation==HOFSApproximationEnum){
     406                        //Do nothing condensatino already done for Stokes part
     407                }
     408                else{
     409                        /*Condense if requested*/
     410                        if(this->element_type==MINIcondensedEnum){
     411                                int indices[3]={18,19,20};
     412                                Ke->StaticCondensation(3,&indices[0]);
     413                        }
     414                        else if(this->element_type==P1bubblecondensedEnum){
     415                                int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     416                                int offset = 0;
     417                                for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     418                                int* indices=xNew<int>(size);
     419                                for(int i=0;i<size;i++) indices[i] = offset+i;
     420                                Ke->StaticCondensation(size,indices);
     421                                xDelete<int>(indices);
     422                        }
    416423                }
    417424
    … …  
    565572                /*StaticCondensation if requested*/
    566573                if(this->element_type==MINIcondensedEnum){
    567                         int indices[3]={18,19,20};
    568 
    569                         this->element_type=MINIEnum;
    570                         ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
    571                         this->element_type=MINIcondensedEnum;
    572 
    573                         pe->StaticCondensation(Ke,3,&indices[0]);
    574                         delete Ke;
     574                        int approximation;
     575                        inputs->GetInputValue(&approximation,ApproximationEnum);
     576                        if(approximation==HOFSApproximationEnum){
     577                                //Do nothing, condensation already done in PVectorCoupling
     578                        }
     579                        else{
     580                                int indices[3]={18,19,20};
     581
     582                                this->element_type=MINIEnum;
     583                                ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
     584                                this->element_type=MINIcondensedEnum;
     585
     586                                pe->StaticCondensation(Ke,3,&indices[0]);
     587                                delete Ke;
     588                        }
     589
    575590                }
    576591                else if(this->element_type==P1bubblecondensedEnum){
    … …  
    66246639        int pnumnodes = this->NumberofNodesPressure();
    66256640
    6626         De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
     6641        De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    66276642
    66286643        for(int i=0;i<vnumnodes;i++){
    … …  
    70687083
    70697084        /*Constants*/
    7070         const int numnodes  = 2 *NUMVERTICES;
    7071         const int numdofp     = NDOF2 *NUMVERTICES;
    7072         const int numdofs     = NDOF4 *NUMVERTICES;
    7073         const int numdoftotal = (NDOF2+NDOF4) *NUMVERTICES;
     7085        const int numnodes       = 2 *NUMVERTICES;
     7086        const int numdofp        = NDOF2 *NUMVERTICES;
     7087        const int numdofs        = NDOF4 * 6 + NDOF3;
     7088        const int numdoftotal    = (NDOF2+NDOF4) *NUMVERTICES + NDOF3;
    70747089
    70757090        /*Intermediaries*/
    7076         Node     *node_list[numnodes];
    7077         int       cs_list[numnodes];
    7078         int       i,j;
    7079         int       init;
    7080 
    7081         /*Prepare node list*/
    7082         for(i=0;i<NUMVERTICES;i++){
    7083                 node_list[i+0*NUMVERTICES] = this->nodes[i];
    7084                 node_list[i+1*NUMVERTICES] = this->nodes[i];
    7085                 cs_list[i+0*NUMVERTICES] = XYEnum;
    7086                 cs_list[i+1*NUMVERTICES] = XYZEnum;
    7087         }
     7091        int       i,j,init;
    70887092
    70897093        /*Some parameters needed*/
    … …  
    70927096        /*compute all stiffness matrices for this element*/
    70937097        ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
    7094         ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
     7098        ElementMatrix* Ke2=new ElementMatrix(this->nodes,2*NUMVERTICES+1,this->parameters,FSvelocityEnum);
    70957099        ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
    70967100        delete Ke1;
    70977101        delete Ke2;
    70987102       
     7103        /*Compute HO Matrix with P1 element type\n");*/
    70997104        this->element_type=P1Enum;
    7100         Ke1=CreateKMatrixDiagnosticHO(); 
     7105        Ke1=CreateKMatrixDiagnosticHO();
    71017106        this->element_type=init;
     7107        /*Compute FS Matrix and condense it \n");*/
    71027108        Ke2=CreateKMatrixDiagnosticFS();
     7109        int indices[3]={18,19,20};
     7110        Ke2->StaticCondensation(3,&indices[0]);
    71037111
    71047112        for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
    7105                 Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF4*j+0];
    7106                 Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF4*j+1];
     7113                Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF3*j+0];
     7114                Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF3*j+1];
    71077115        }
    71087116        for(i=0;i<numdofp;i++) for(j=0;j<NUMVERTICES;j++){
    7109                 Ke->values[i*numdoftotal+numdofp+NDOF4*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0];
    7110                 Ke->values[i*numdoftotal+numdofp+NDOF4*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1];
    7111         }
    7112 
    7113         /*Transform Coordinate System*/
    7114         TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
     7117                Ke->values[i*numdoftotal+numdofp+NDOF3*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0];
     7118                Ke->values[i*numdoftotal+numdofp+NDOF3*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1];
     7119        }
     7120
     7121        /*Transform Coordinate System*/ //Do not transform, already sone in the matrices
     7122        //TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
    71157123
    71167124        /*clean-up and return*/
    … …  
    71197127        return Ke;
    71207128}
    7121 /*}}}*/
     7129//*}}}*/
    71227130/*FUNCTION Penta::CreateKMatrixDiagnosticHoriz {{{*/
    71237131ElementMatrix* Penta::CreateKMatrixDiagnosticHoriz(void){
    … …  
    76577665
    76587666        /*compute all stiffness matrices for this element*/
     7667        ElementMatrix* Ke1=CreateKMatrixDiagnosticFS();
    76597668        int init = this->element_type;
    7660         this->element_type=P1Enum;
    7661         ElementMatrix* Ke1=CreateKMatrixDiagnosticHO();
     7669        this->element_type=P1Enum; //P1 needed for HO
     7670        ElementMatrix* Ke2=CreateKMatrixDiagnosticHO();
    76627671        this->element_type=init;
    7663         ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
    76647672        ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
    76657673        ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
    … …  
    82308238        inputs->GetInputValue(&approximation,ApproximationEnum);
    82318239        if(approximation!=HOFSApproximationEnum) return NULL;
    8232         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
     8240        int vnumnodes = this->NumberofNodesVelocity();
     8241        int pnumnodes = this->NumberofNodesPressure();
     8242
     8243        /*Prepare coordinate system list*/
     8244        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     8245        for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
     8246        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     8247
     8248        /*Initialize Element matrix and vectors*/
     8249        ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    82338250
    82348251        /*Retrieve all inputs and parameters*/
    82358252        GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    82368253        this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
    8237         Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    8238         Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    8239         Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
     8254        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     8255        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     8256        Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    82408257        Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
    82418258
    … …  
    82568273
    82578274                for(i=0;i<NUMVERTICES;i++){
    8258                         pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
    8259                         pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    8260                         pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    8261                         pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
     8275                        pe->values[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
     8276                        pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
     8277                        pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
     8278                        pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    82628279                }
    82638280        }
    82648281
    82658282        /*Transform coordinate system*/
    8266         TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
     8283        TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    82678284
    82688285        /*Clean up and return*/
     8286        xDelete<int>(cs_list);
    82698287        delete gauss;
    82708288        return pe;
    … …  
    84078425
    84088426        /*compute all load vectors for this element*/
     8427        int init = this->element_type;
     8428        this->element_type=P1Enum;
    84098429        ElementVector* pe1=CreatePVectorDiagnosticHO();
     8430        this->element_type=init;
    84108431        ElementVector* pe2=CreatePVectorDiagnosticFS();
     8432        int indices[3]={18,19,20};
     8433        this->element_type=MINIcondensedEnum;
     8434        ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
     8435        this->element_type=init;
     8436        pe2->StaticCondensation(Ke,3,&indices[0]);
     8437        delete Ke;
    84118438        ElementVector* pe3=CreatePVectorCouplingHOFS();
    8412         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     8439        ElementVector* pe =new ElementVector(pe1,pe2);
    84138440
    84148441        /*clean-up and return*/

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp

    r15696 r15744  
    100100                        for(int i=0;i<iomodel->numberofelements;i++){
    101101                                if(iomodel->my_elements[i]){
    102                                         nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,DiagnosticHorizAnalysisEnum,FSvelocityEnum));
     102                                        Node* node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,DiagnosticHorizAnalysisEnum,FSvelocityEnum);
     103                                        node->Deactivate();
     104                                        nodes->AddObject(node);
    103105                                }
    104106                        }
Note: See TracChangeset for help on using the changeset viewer.