Changeset 5938


Ignore:
Timestamp:
09/22/10 09:15:54 (15 years ago)
Author:
Mathieu Morlighem
Message:

more Element Matrix and Vector

Location:
issm/trunk/src/c/objects/Loads
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Loads/Riftfront.cpp

    r5911 r5938  
    1919#include "../objects.h"
    2020/*}}}*/
     21
     22/*Element macros*/
     23#define NUMVERTICES 2
     24#define NDOF1 1
     25#define NDOF2 2
     26#define NDOF3 3
     27#define NDOF4 4
    2128
    2229/*Riftfront constructors and destructor*/
    … …  
    375382}
    376383/*}}}*/
     384/*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/
     385void  Riftfront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
     386
     387        /*Retrieve parameters: */
     388        ElementMatrix* Ke=NULL;
     389        int analysis_type;
     390        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     391
     392        switch(analysis_type){
     393                case DiagnosticHorizAnalysisEnum:
     394                        Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax);
     395                        break;
     396                case AdjointHorizAnalysisEnum:
     397                        Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax);
     398                        break;
     399                default:
     400                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     401        }
     402
     403        /*Add to global Vector*/
     404        if(Ke){
     405                Ke->AddToGlobal(Kgg,Kff,Kfs);
     406                delete Ke;
     407        }
     408}
     409/*}}}1*/
     410/*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/
     411void  Riftfront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
     412
     413        /*Retrieve parameters: */
     414        ElementVector* pe=NULL;
     415        int analysis_type;
     416        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     417
     418        switch(analysis_type){
     419                case DiagnosticHorizAnalysisEnum:
     420                        pe=PenaltyCreatePVectorDiagnosticHoriz(kmax);
     421                        break;
     422                case AdjointHorizAnalysisEnum:
     423                        /*No penalty applied on load vector*/
     424                        break;
     425                default:
     426                        ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
     427        }
     428
     429        /*Add to global Vector*/
     430        if(pe){
     431                pe->AddToGlobal(pg,pf);
     432                delete pe;
     433        }
     434}
     435/*}}}1*/
    377436/*FUNCTION Riftfront::CreateKMatrix {{{1*/
    378437void  Riftfront::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){
    … …  
    384443void  Riftfront::CreatePVector(Vec pg,Vec pf){
    385444        /*do nothing: */
    386 }
    387 /*}}}1*/
    388 /*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/
    389 void  Riftfront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){
    390 
    391         int         i;
    392         int         j;
    393         const int   numgrids            = MAX_RIFTFRONT_GRIDS;
     445        return;
     446}
     447/*}}}1*/
     448/*FUNCTION Riftfront::InAnalysis{{{1*/
     449bool Riftfront::InAnalysis(int in_analysis_type){
     450        if (in_analysis_type==this->analysis_type) return true;
     451        else return false;
     452}
     453/*}}}*/
     454
     455/*Riftfront numerics*/
     456/*FUNCTION Riftfront::PenaltyCreateKMatrixDiagnosticHoriz {{{1*/
     457ElementMatrix* Riftfront::PenaltyCreateKMatrixDiagnosticHoriz(double kmax){
     458
     459        const int   numdof = NDOF2*NUMVERTICES;
     460        int         i,j;
    394461        int         dofs[1]             = {0};
    395462        double      Ke_gg[4][4];
    396         const int   numdof              = 2 *numgrids;
    397         int*        doflist=NULL;
    398463        double      thickness;
    399464        double      h[2];
    … …  
    413478        /*enum of element? */
    414479        if(elements[0]->Enum()!=TriaEnum)ISSMERROR(" only Tria element allowed for Riftfront load!");
    415 
    416         /*recover elements on both side of rift: */
    417480        tria1=(Tria*)elements[0];
    418481        tria2=(Tria*)elements[1];
    419482
    420        
    421         /* Get node coordinates and dof list: */
    422         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    423 
    424         /* Set Ke_gg to 0: */
    425         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
     483        /*Initialize Element Matrix*/
     484        if(!this->active) return NULL;
     485        ElementMatrix* Ke=NewElementMatrix(nodes,NUMVERTICES,this->parameters);
    426486
    427487        /*Get some parameters: */
    428488        this->parameters->FindParam(&penalty_offset,PenaltyOffsetEnum);
    429489        this->inputs->GetParameterValue(&friction,FrictionEnum);
    430 
    431         if(this->active){
    432        
    433                 /*There is contact, we need to constrain the normal velocities (zero penetration), and the
    434                  *contact slip friction. */
    435                  
    436                 /*Recover thickness: */
    437                 tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum);
    438                 tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum);
    439 
    440                 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
    441                 thickness=h[0];
    442 
    443                 /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
    444                 //First line:
    445                 Ke_gg[0][0]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
    446                 Ke_gg[0][1]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    447                 Ke_gg[0][2]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
    448                 Ke_gg[0][3]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    449                 //Second line:
    450                 Ke_gg[1][0]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    451                 Ke_gg[1][1]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
    452                 Ke_gg[1][2]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    453                 Ke_gg[1][3]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
    454                 //Third line:
    455                 Ke_gg[2][0]+=-pow(normal[0],2)*kmax*pow(10,penalty_offset);
    456                 Ke_gg[2][1]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    457                 Ke_gg[2][2]+=pow(normal[0],2)*kmax*pow(10,penalty_offset);
    458                 Ke_gg[2][3]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    459                 //Fourth line:
    460                 Ke_gg[3][0]+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    461                 Ke_gg[3][1]+=-pow(normal[1],2)*kmax*pow(10,penalty_offset);
    462                 Ke_gg[3][2]+=normal[0]*normal[1]*kmax*pow(10,penalty_offset);
    463                 Ke_gg[3][3]+=pow(normal[1],2)*kmax*pow(10,penalty_offset);
    464 
    465                 /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
    466                
    467                 //First line:
    468                 Ke_gg[0][0]+=pow(normal[1],2)*thickness*length*friction;
    469                 Ke_gg[0][1]+=-normal[0]*normal[1]*thickness*length*friction;
    470                 Ke_gg[0][2]+=-pow(normal[1],2)*thickness*length*friction;
    471                 Ke_gg[0][3]+=normal[0]*normal[1]*thickness*length*friction;
    472                 //Second line:
    473                 Ke_gg[1][0]+=-normal[0]*normal[1]*thickness*length*friction;
    474                 Ke_gg[1][1]+=pow(normal[0],2)*thickness*length*friction;
    475                 Ke_gg[1][2]+=normal[0]*normal[1]*thickness*length*friction;
    476                 Ke_gg[1][3]+=-pow(normal[0],2)*thickness*length*friction;
    477                 //Third line:
    478                 Ke_gg[2][0]+=-pow(normal[1],2)*thickness*length*friction;
    479                 Ke_gg[2][1]+=normal[0]*normal[1]*thickness*length*friction;
    480                 Ke_gg[2][2]+=pow(normal[1],2)*thickness*length*friction;
    481                 Ke_gg[2][3]+=-normal[0]*normal[1]*thickness*length*friction;
    482                 //Fourth line:
    483                 Ke_gg[3][0]+=normal[0]*normal[1]*thickness*length*friction;
    484                 Ke_gg[3][1]+=-pow(normal[0],2)*thickness*length*friction;
    485                 Ke_gg[3][2]+=-normal[0]*normal[1]*thickness*length*friction;
    486                 Ke_gg[3][3]+=pow(normal[0],2)*thickness*length*friction;
    487 
    488                 /*Add Ke_gg to global matrix Kgg: */
    489                 MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    490         }
    491         else{
    492                 /*the grids on both sides of the rift do not penetrate.  PenaltyCreatePVector will
    493                  *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness,
    494                  there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */
    495         }
    496         /*Free ressources:*/
    497         xfree((void**)&doflist);
    498 
    499 }
    500 /*}}}1*/
    501 /*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/
    502 void  Riftfront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){
    503 
    504         int         i                     ,j;
    505         const int   numgrids            = MAX_RIFTFRONT_GRIDS;
    506         double      pe_g[4]={0.0};
    507         const int   numdof              = 2 *numgrids;
    508         int*        doflist=NULL;
    509 
     490        tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum);
     491        tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum);
     492        if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
     493        thickness=h[0];
     494
     495        /*There is contact, we need to constrain the normal velocities (zero penetration), and the
     496         *contact slip friction. */
     497
     498        /*From Peter Wriggers book (Computational Contact Mechanics, p191): */
     499        Ke->values[0*numdof+0]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
     500        Ke->values[0*numdof+1]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     501        Ke->values[0*numdof+2]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
     502        Ke->values[0*numdof+3]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     503
     504        Ke->values[1*numdof+0]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     505        Ke->values[1*numdof+1]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
     506        Ke->values[1*numdof+2]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     507        Ke->values[1*numdof+3]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
     508
     509        Ke->values[2*numdof+0]+= -pow(normal[0],2)*kmax*pow(10,penalty_offset);
     510        Ke->values[2*numdof+1]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     511        Ke->values[2*numdof+2]+= +pow(normal[0],2)*kmax*pow(10,penalty_offset);
     512        Ke->values[2*numdof+3]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     513
     514        Ke->values[3*numdof+0]+= -normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     515        Ke->values[3*numdof+1]+= -pow(normal[1],2)*kmax*pow(10,penalty_offset);
     516        Ke->values[3*numdof+2]+= +normal[0]*normal[1]*kmax*pow(10,penalty_offset);
     517        Ke->values[3*numdof+3]+= +pow(normal[1],2)*kmax*pow(10,penalty_offset);
     518
     519        /*Now take care of the friction: of type sigma=frictiontangent_velocity2-tangent_velocity1)*/
     520
     521        Ke->values[0*numdof+0]+= +pow(normal[1],2)*thickness*length*friction;
     522        Ke->values[0*numdof+1]+= -normal[0]*normal[1]*thickness*length*friction;
     523        Ke->values[0*numdof+2]+= -pow(normal[1],2)*thickness*length*friction;
     524        Ke->values[0*numdof+3]+= +normal[0]*normal[1]*thickness*length*friction;
     525
     526        Ke->values[1*numdof+0]+= -normal[0]*normal[1]*thickness*length*friction;
     527        Ke->values[1*numdof+1]+= +pow(normal[0],2)*thickness*length*friction;
     528        Ke->values[1*numdof+2]+= +normal[0]*normal[1]*thickness*length*friction;
     529        Ke->values[1*numdof+3]+= -pow(normal[0],2)*thickness*length*friction;
     530
     531        Ke->values[2*numdof+0]+= -pow(normal[1],2)*thickness*length*friction;
     532        Ke->values[2*numdof+1]+= +normal[0]*normal[1]*thickness*length*friction;
     533        Ke->values[2*numdof+2]+= +pow(normal[1],2)*thickness*length*friction;
     534        Ke->values[2*numdof+3]+= -normal[0]*normal[1]*thickness*length*friction;
     535
     536        Ke->values[3*numdof+0]+= +normal[0]*normal[1]*thickness*length*friction;
     537        Ke->values[3*numdof+1]+= -pow(normal[0],2)*thickness*length*friction;
     538        Ke->values[3*numdof+2]+= -normal[0]*normal[1]*thickness*length*friction;
     539        Ke->values[3*numdof+3]+= +pow(normal[0],2)*thickness*length*friction;
     540
     541        /*Clean up and return*/
     542        return Ke;
     543}
     544/*}}}1*/
     545/*FUNCTION Riftfront::PenaltyCreatePVectorDiagnosticHoriz {{{1*/
     546ElementVector* Riftfront::PenaltyCreatePVectorDiagnosticHoriz(double kmax){
     547
     548        const int   numdof = NDOF2*NUMVERTICES;
     549        int         i,j;
    510550        double      rho_ice;
    511551        double      rho_water;
    … …  
    523563        bool        shelf;
    524564
    525 
    526565        /*Objects: */
    527566        Element   **elements            = NULL;
    … …  
    531570        Matpar     *matpar              = NULL;
    532571
    533        
    534572        /*Recover hook objects: */
    535573        elements=(Element**)helements->deliverp();
    … …  
    539577        /*enum of element? */
    540578        if(elements[0]->Enum()!=TriaEnum)ISSMERROR(" only Tria element allowed for Riftfront load!");
    541 
    542         /*recover elements on both side of rift: */
    543579        tria1=(Tria*)elements[0];
    544580        tria2=(Tria*)elements[1];
    545581
    546         /* Get node coordinates and dof list: */
    547         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     582        /*Initialize Element Matrix*/
     583        if(this->active) return NULL; /*The penalty is active. No loads implied here.*/
     584        ElementVector* pe=NewElementVector(nodes,NUMVERTICES,this->parameters);
    548585
    549586        /*Get some inputs: */
    550587        this->inputs->GetParameterValue(&fill,FillEnum);
    551588        this->inputs->GetParameterValue(&shelf,SegmentOnIceShelfEnum);
    552 
    553         if(!this->active){
    554                 /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
    555                  * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/
    556        
    557                 /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to
    558                  * be the same across the rift.: */
    559 
    560                 rho_ice=matpar->GetRhoIce();
    561                 rho_water=matpar->GetRhoWater();
    562                 gravity=matpar->GetG();
    563 
    564                 /*get thickness: */
    565                 tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum);
    566                 tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum);
    567 
    568                 if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
    569                 thickness=h[0];
    570 
    571                 tria1->GetParameterValue(&b[0],nodes[0],BedEnum);
    572                 tria2->GetParameterValue(&b[1],nodes[1],BedEnum);
    573 
    574                 if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts");
    575                 bed=b[0];
    576 
    577                 /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
    578                 if(fill==WaterEnum){
    579                         if(shelf){
    580                                 /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
    581                                 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2  - rho_water*gravity*pow(bed,(double)2)/(double)2;
    582                         }
    583                         else{
    584                                 //We are on an icesheet, we assume the water column fills the entire front: */
    585                                 pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2  - rho_water*gravity*pow(thickness,(double)2)/(double)2;
    586                         }
    587                 }
    588                 else if(fill==AirEnum){
    589                         pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2;   //icefront on an ice sheet, pressure imbalance ice vs air.
    590                 }
    591                 else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
    592                         pressure=0;
    593                 }
    594                 else if(fill==MelangeEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
    595                        
    596                         if(!shelf) ISSMERROR("%s%i%s","fill type ",fill," not supported on ice sheets yet.");
    597 
    598                         pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2;
    599                         pressure_air=0;
    600                         pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2;
    601                         pressure_water=1.0/2.0*rho_water*gravity*  ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
    602 
    603                         pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
     589        rho_ice=matpar->GetRhoIce();
     590        rho_water=matpar->GetRhoWater();
     591        gravity=matpar->GetG();
     592        tria1->GetParameterValue(&h[0],nodes[0],ThicknessEnum);
     593        tria2->GetParameterValue(&h[1],nodes[1],ThicknessEnum);
     594        if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts");
     595        thickness=h[0];
     596        tria1->GetParameterValue(&b[0],nodes[0],BedEnum);
     597        tria2->GetParameterValue(&b[1],nodes[1],BedEnum);
     598        if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts");
     599        bed=b[0];
     600
     601        /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics,
     602         * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/
     603
     604        /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to
     605         * be the same across the rift.: */
     606
     607        /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */
     608        if(fill==WaterEnum){
     609                if(shelf){
     610                        /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */
     611                        pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2  - rho_water*gravity*pow(bed,(double)2)/(double)2;
    604612                }
    605613                else{
    606                         ISSMERROR("%s%i%s","fill type ",fill," not supported yet.");
     614                        //We are on an icesheet, we assume the water column fills the entire front: */
     615                        pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2  - rho_water*gravity*pow(thickness,(double)2)/(double)2;
    607616                }
    608 
    609                 /*Ok, add contribution to first grid, along the normal i==0: */
    610                 for (j=0;j<2;j++){
    611                         pe_g[j]+=pressure*normal[j]*length;
    612                 }
    613        
    614                 /*Add contribution to second grid, along the opposite normal: i==1 */
    615                 for (j=0;j<2;j++){
    616                         pe_g[2+j]+= -pressure*normal[j]*length;
    617                 }       
    618                 /*Add pe_g to global vector pg; */
    619                 VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    620 
     617        }
     618        else if(fill==AirEnum){
     619                pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2;   //icefront on an ice sheet, pressure imbalance ice vs air.
     620        }
     621        else if(fill==IceEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
     622                pressure=0;
     623        }
     624        else if(fill==MelangeEnum){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice)
     625
     626                if(!shelf) ISSMERROR("%s%i%s","fill type ",fill," not supported on ice sheets yet.");
     627
     628                pressure_litho=rho_ice*gravity*pow(thickness,(double)2)/(double)2;
     629                pressure_air=0;
     630                pressure_melange=rho_ice*gravity*pow(fraction*thickness,(double)2)/(double)2;
     631                pressure_water=1.0/2.0*rho_water*gravity*  ( pow(bed,2.0)-pow(rho_ice/rho_water*fraction*thickness,2.0) );
     632
     633                pressure=pressure_litho-pressure_air-pressure_melange-pressure_water;
    621634        }
    622635        else{
    623                 /*The penalty is active. No loads implied here.*/
    624         }
    625        
    626         /*Free ressources:*/
    627         xfree((void**)&doflist);
    628 
    629 }
    630 /*}}}1*/
    631 /*FUNCTION Riftfront::InAnalysis{{{1*/
    632 bool Riftfront::InAnalysis(int in_analysis_type){
    633         if (in_analysis_type==this->analysis_type) return true;
    634         else return false;
    635 }
    636 /*}}}*/
    637 
    638 /*Riftfront numerics*/
     636                ISSMERROR("%s%i%s","fill type ",fill," not supported yet.");
     637        }
     638
     639        /*Ok, add contribution to first grid, along the normal i==0: */
     640        for (j=0;j<2;j++){
     641                pe->values[j]+=pressure*normal[j]*length;
     642        }
     643
     644        /*Add contribution to second grid, along the opposite normal: i==1 */
     645        for (j=0;j<2;j++){
     646                pe->values[2+j]+= -pressure*normal[j]*length;
     647        }       
     648
     649        /*Clean up and return*/
     650        return pe;
     651}
     652/*}}}1*/
    639653/*FUNCTION Riftfront::Constrain {{{1*/
    640654#define _ZIGZAGCOUNTER_
    … …  
    731745}
    732746/*}}}1*/
    733 /*FUNCTION Riftfront::GetDofList {{{1*/
    734 void  Riftfront::GetDofList(int** pdoflist,int approximation,int setenum){
    735 
    736         int i,j;
    737         int numberofdofs=0;
    738         int count=0;
    739 
    740         /*output: */
    741         int* doflist=NULL;
    742 
    743         /*pointers: */
    744         Node**   nodes=NULL;
    745        
    746         /*recover pointers: */
    747         nodes=(Node**)hnodes->deliverp();
    748        
    749         /*Some checks for debugging*/
    750         ISSMASSERT(nodes);
    751                
    752         /*Figure out size of doflist: */
    753         for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
    754                 numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    755         }
    756 
    757         /*Allocate: */
    758         doflist=(int*)xmalloc(numberofdofs*sizeof(int));
    759 
    760         /*Populate: */
    761         count=0;
    762         for(i=0;i<MAX_RIFTFRONT_GRIDS;i++){
    763                 nodes[i]->GetDofList(doflist+count,approximation,setenum);
    764                 count+=nodes[i]->GetNumberOfDofs(approximation,setenum);
    765         }
    766 
    767         /*Assign output pointers:*/
    768         *pdoflist=doflist;
    769 }
    770 /*}}}*/
    771747/*FUNCTION Riftfront::IsFrozen{{{1*/
    772748bool   Riftfront::IsFrozen(void){

  • issm/trunk/src/c/objects/Loads/Riftfront.h

    r5772 r5938  
    1313class Inputs;
    1414class IoModel;
    15 
    16 #define MAX_RIFTFRONT_GRIDS 2 //max number of grids on a rift flank, only 2 because 2d for now.
    17 #define RIFTFRONTSTRING 20 //max string length
    1815/*}}}*/
    1916
    … …  
    8178                /*}}}*/
    8279                /*Riftfront specific routines: {{{1*/
    83                 void  GetDofList(int** pdoflist,int approximation_enum,int setenum);
    8480                bool  PreStable();
     81                ElementMatrix* PenaltyCreateKMatrixDiagnosticHoriz(double kmax);
     82                ElementVector* PenaltyCreatePVectorDiagnosticHoriz(double kmax);
    8583                void  SetPreStable();
    8684                int   PreConstrain(int* punstable);
Note: See TracChangeset for help on using the changeset viewer.