Index: /issm/trunk/src/c/objects/cielo/PenpairLoad.c =================================================================== --- /issm/trunk/src/c/objects/cielo/PenpairLoad.c (revision 215) +++ /issm/trunk/src/c/objects/cielo/PenpairLoad.c (revision 215) @@ -0,0 +1,1505 @@ + +/* + PenpairLoad.c: +*/ + +/* PenpairLoad: this object deals with the connection between a 2d and a + * 3d formulation, for a pair of grids. + * Given a 2d grid A, with a formulation involving for example dofs + * 1 and 2 (ex: horizontal velocity), and another grid B belonging to the + * 3d mesh, we want to constrain one of the dofs, for example 1, so that + * values for this dof are the same on grid A and on grid B. + * This is feasible by applying a penalty stiffness matrix + * for this pair of grids, for the right dof we are dealing with.*/ + +#include +#include +#include +#include + +/* environment: */ +#include "../../include/cielo_types.h" +#include "../../include/cielo_macros.h" +#include "../../include/cielo_externals.h" + +/* hosting environment: */ +#include "../../include/matlab_includes.h" + +/* object handlers: */ +#include "../../DataSet/DataSet.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +/* shared utilities (includes alloc.h memory management): */ +#include "../../shared/shared.h" + +/* "virtual function" access: */ +#include "./../VirtualFunctions.h" + +/* object itself: */ +#include "./PenpairLoad.h" + +/*Objects: */ +#include "../objects.h" + +#define NDOF2 2 + +/*-------------------------------------------------- + NewPenpairLoad + --------------------------------------------------*/ + +PenpairLoad* NewPenpairLoad( void* vppenpair) +{ + + Penpair* penpair = (Penpair*)vppenpair; + PenpairLoad* this = NULL; + int i,j; + + this=(PenpairLoad*)xcalloc(1,sizeof(PenpairLoad)); /* guaranteed NULL, 0 initialisation */ + + /* include the "base class" in the "derived" object: */ + memcpy( &this->penpair, penpair, sizeof(Penpair)); + + /* since zero could be a valid index, set all DataSet index + "pointers" to UNDEF: */ + for ( i=0; i<2; i++) this->internalgrid_indices[i] = UNDEF; + for ( i=0; i<2; i++) this->internalgrid_indicesb[i] = UNDEF; + for ( i=0; i<2; i++) this->pg[i]=NULL; + + for ( i=0; i<2; i++) this->internalelement_indices[i] = UNDEF; + for ( i=0; i<2; i++) this->elements[i]=NULL; + + this->prestable=0; //set of constraints is not stable in the beginning + this->active=0; //no constraints active at the beginning. + this->counter=0; + + return this; +} + + +/*-------------------------------------------------- + DeletePenpairLoad + --------------------------------------------------*/ + +void DeletePenpairLoad(void* *vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + PenpairLoad* *pthis = (PenpairLoad**)vpthis; + PenpairLoad* this = *pthis; + + xfree(vpthis); +} + + +/*-------------------------------------------------- + PenpairLoadEcho + --------------------------------------------------*/ + +void PenpairLoadEcho(void* vpthis){ + int i; + + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = (PenpairLoad*)vpthis; + + if (this) { + + printf("\nPenpairLoad echo:\n"); + printf("-------------------\n"); + printf(" base Penpair object:\n"); + PenpairEcho( &this->penpair); + printf("\n"); + printf(" penalty pair points to following grids (pointer-based ops):\n"); + printf(" local grid number corresponding grid id\n"); + printf(" ----------------- ---------------------\n"); + for ( i=0; i<2; i++) { + if ( this->pg[i]) + printf(" %8d %8d\n", i, InternalGridGetID( this->pg[i])); + } + printf(" penalty pair points to following elements (pointer-based ops):\n"); + printf(" local element number corresponding element id\n"); + printf(" ----------------- ---------------------\n"); + for ( i=0; i<2; i++) { + if ( this->elements[i]){ + TriaElement* triaelement=this->elements[i]; + printf(" %8d %8d\n", i, TriaElementGetID(triaelement)); + } + } + printf(" active %i \n",this->active); + printf(" counter %i \n",this->counter); + printf(" prestable %i \n",this->prestable); + + } + else{ + printf(" null pointer input for PenpairLoad object; no echo\n"); + } +} + +/*-------------------------------------------------- + PenpairLoadSizeOfInDoubles + --------------------------------------------------*/ + +int PenpairLoadSizeOfInDoubles(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + PenpairLoad* this = (PenpairLoad*)vpthis; + int struct_doubles=0; + + struct_doubles = + (int)ceil((float)sizeof(*this)/(float)sizeof(double)); + + #ifdef _DEBUG_ + printf("\n PenpairLoadSizeOfInDoubles diagnostics:\n"); + printf(" returning %d doubles\n", struct_doubles); + #endif + + return struct_doubles; +} + + +/*-------------------------------------------------- + PenpairLoadShrinkWrap + --------------------------------------------------*/ + +void PenpairLoadShrinkWrap(void* to, void* vpthis) +{ + /* no "extended" data means a simple copy will do: */ + + memcpy( to, vpthis, sizeof(double)*PenpairLoadSizeOfInDoubles(vpthis)); +} + + +/*-------------------------------------------------- + PenpairLoadUnwrap + --------------------------------------------------*/ + +int PenpairLoadUnwrap(void* vpthis) +{ + /* no extended data means no internal pointers to reestablish */ + + return 1; +} + +/*-------------------------------------------------- + PenpairLoadMarshall + --------------------------------------------------*/ + +void PenpairLoadMarshall( char* *pbuffer, void* vpthis, int* size) +{ + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = (PenpairLoad*)vpthis; + + char* buffer=*pbuffer; + + int i, j, size_penpair=0; + int ig; + + int num_ints=0, num_doubles=0; + + + /* first marshall the embedded Penpair object: ... */ + + PenpairMarshall( &buffer, &this->penpair, &size_penpair); + + + /* ... then marshall the PenpairLoad member data: */ + + for ( i=0; i<2; i++) { + cellmar( &this->internalgrid_indices[i], &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->internalgrid_indicesb[i], &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->pg[i], &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->internalelement_indices[i], &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->elements[i], &buffer, INTEGER_FLAG); num_ints++; + + } + cellmar( &this->active, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->counter, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->prestable, &buffer, INTEGER_FLAG); num_ints++; + + *size = size_penpair + num_ints*sizeof(int) + num_doubles*sizeof(double); + *pbuffer = buffer; + + return; +} + + +/*-------------------------------------------------- + PenpairLoadDemarshall + --------------------------------------------------*/ + +void* PenpairLoadDemarshall( DataSet* ds, char* *pbuffer, int machine_flag) +{ + int i, j, ig; + Penpair* ppenpair=NULL; /* for embedded object */ + + char* buffer=*pbuffer; + PenpairLoad* this=NULL; + + + /* first demarshall the embedded Penpair object: ... */ + + ppenpair = PenpairDemarshall( NULL, &buffer, machine_flag); + /* create the derived object: */ + this = NewPenpairLoad( ppenpair); + /* and free penpair now that it's embedded in the derived object: */ + DeletePenpair( (void**)&ppenpair); + + + /* ... then demarshall the PenpairLoad member data: */ + + for ( i=0; i<2; i++) { + celldem( &this->internalgrid_indices[i], &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->internalgrid_indicesb[i], &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->pg[i], &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->internalgrid_indices[i], &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->internalelement_indices[i], &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->elements[i], &buffer, INTEGER_FLAG, machine_flag); + + } + celldem( &this->active, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->counter, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->prestable, &buffer, INTEGER_FLAG, machine_flag); + + /*Set pointers to NULL :*/ + for (i=0;i<2;i++){ + this->pg[i]=NULL; + this->elements[i]=NULL; + } + + if (ds) PenpairLoadAddToDataSet( ds, this); + *pbuffer=buffer; + + return (void*)this; +} + +/*-------------------------------------------------- + PenpairLoadGetID + --------------------------------------------------*/ + +int PenpairLoadGetID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = (PenpairLoad*)vpthis; + + return PenpairGetID( &this->penpair); +} + +/*-------------------------------------------------- + PenpairLoadSetVirtualFunctions + --------------------------------------------------*/ + +int PenpairLoadSetVirtualFunctions( void* vp) +{ + + /* VirtualFunctions implemented as a "friend" class: */ + + /* void* for polymorphic function compatibility */ + VirtualFunctions* vf = (VirtualFunctions*)vp; + + if ( VirtualFunctionsInit(vf) ) { + + /* set PenpairLoad-specific virtual function pointers: */ + + /* general: */ + + /* configuration: */ + vf->Configure = &PenpairLoadConfigure; + + /* penalty stiffnesses: */ + vf->CreateKMatrix = NULL; + vf->PenaltyCreateKMatrix = &PenpairLoadPenaltyCreateKMatrix; + vf->PenaltyCreatePVector = &PenpairLoadPenaltyCreatePVector; + vf->PenaltyConstrain = &PenpairLoadPenaltyConstrain; + vf->PenaltyPreConstrain = &PenpairLoadPenaltyPreConstrain; + vf->PenaltyPotentialUnstableConstraint = &PenpairLoadPotentialUnstableConstraint; + vf->Penetration = &PenpairLoadPenetration; + + return 1; + } + else + return 0; + +} + +/*-------------------------------------------------- + PenpairLoadAddToDataSet + --------------------------------------------------*/ + +int PenpairLoadAddToDataSet( DataSet* dataset, PenpairLoad* this) +{ + /* + + Provided as an alternate interface to DataSetAddObject. + + In the c++ world, the dataset could just add an object to its list + using a dataset member function and a pointer to the object's base + class. + + In the c world, however, only the object itself is capable of telling + the dataset its type, hence the structure here; a call to a generic + DataSetAddObject function via an object "member function." + + */ + + DataSetAddObject( dataset, PenpairLoadEnum(), (void*)this, + &DeletePenpairLoad, &PenpairLoadEcho, &PenpairLoadSizeOfInDoubles, + &PenpairLoadShrinkWrap, &PenpairLoadUnwrap, + &PenpairLoadMarshall, + &PenpairLoadGetID, &PenpairLoadSetVirtualFunctions); + + return 1; +} + +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadConfigure" + +/*-------------------------------------------------- + PenpairLoadConfigure + --------------------------------------------------*/ + +int PenpairLoadConfigure( void* vpthis, int num_datasets, ...) { + + /* set up data connections to all necessary, related objects: */ + + PenpairLoad* this = (PenpairLoad*)vpthis; + + int i,j, found, foundgrid, noerr=1, first=1; + + /* expected DataSet input: */ + + DataSet* pgpdt=NULL; + DataSet* pgpdtb=NULL; + DataSet* est=NULL; + + /* variable argument list handling: */ + + va_list args; + DataSet* *arglist=NULL; + + #ifdef _DEBUG_ + printf("virtual function: PenpairLoadConfigure\n"); + #endif + + + + if (this) { + + + if (!num_datasets) { + + printf(" %s error:\n", __FUNCT__ ); + printf(" null input for number of datasets;\n"); + printf(" could not establish grid connections.\n"); + return 0; + } + else { + + /* read, save dataset pointers for data connection operations: */ + arglist = xmalloc( num_datasets*sizeof(DataSet*)); + va_start( args, num_datasets); + for ( i=0; ipg[i]*/ + for (i=0;i<2;i++){ + this->pg[i]=NULL; + } + for (i=0;i<2;i++){ + this->elements[i]=NULL; + } + + /* grid connections based on DataSet of type "gpdt": */ + + #ifndef _PARALLEL_ + found=0; + for ( i=0; ipenpair); + + for ( i=0; i<2; i++) { + + if ( (current_grid_id=pgrid_id_list[i]) ) { /* assign and test presence */ + + /* if grid pointer information has been established once already + (e.g., by a previous module or operation), grid index information + may already be present and valid. check to see if corresponding + pointer information can be reestablished without an expensive + seek operation: */ + + if ( UNDEF != (int)this->internalgrid_indices[i] ) { + + /* does index still point to a valid grid? (i.e., the one with + the same id): */ + + if ( InternalGridEnum() == DataSetGetEnum( pgpdt, this->internalgrid_indices[i]) + && + current_grid_id == DataSetGetObjectID( pgpdt, this->internalgrid_indices[i]) ) { + + /* great, still valid! since we have the index position in + gpdt, just use it to get, save, the grid pointer: */ + + this->pg[i] = DataSetGetObjectPtr( pgpdt, this->internalgrid_indices[i]); + } + } + + /* if we haven't yet "recovered" the InternalGrid object pointer, + go seek it out: */ + + if ( !this->pg[i] ) { + + noerr = DataSetSeek( (void**)(this->pg+i), /* cast because some compilers recognise as */ + /* ptr to ptr to struct otherwise */ + (this->internalgrid_indices+i) , + pgpdt, current_grid_id, 1, InternalGridEnum()); + + if (!noerr) { + printf(" %s error:\n", __func__); + printf("%s%s%i%s%i%s%s%s\n",__FUNCT__," error message: could not find grid with id of ",current_grid_id," for load ",PenpairLoadGetID(this)," of type \"",PenpairLoadGetName(this),"\""); + } + } + } + } + } + else { + + if (first) { + printf(" %s error:\n", __func__); + first=0; + } + printf("%s%s\n",__FUNCT__, "error message: could not find DataSet of type \"gpdt\", grid object\n"); + printf(" pointers not established for penpair load id = %d\n", + PenpairLoadGetID(this)); + noerr = 0; + } + + + + #else //ifdef _PARALLEL_ + + /*Configure the load when we are running on a parallel cluster: in this case, + we have partitioned the grids and elements across the cluster. Every node has + its own lst, its own bgpdt (which owns the internal grids of the partition) and + a common bgpdtb (b for boundaries of the partition). The boundary grids have to be + dealt with. What happens is that bgpdt might not own all the grids referenced by this + load, bgpdtb might own some of these grids. So we have to do a double configuration + sort of: + */ + + /*First, we need to recover pgpdt and pgpdtb from the arglist: */ + found=0; + for ( i=0; ipenpair); + + for ( i=0; i<2; i++) { + + if ( (current_grid_id=pgrid_id_list[i]) ) { /* assign and test presence */ + + /* if grid pointer information has been established once already + (e.g., by a previous module or operation), grid index information + may already be present and valid. check to see if corresponding + pointer information can be reestablished without an expensive + seek operation: */ + + if ( UNDEF != (int)this->internalgrid_indices[i] ) { + + /* does index still point to a valid grid? (i.e., the one with + the same id): */ + + if ( InternalGridEnum() == DataSetGetEnum( pgpdt, this->internalgrid_indices[i]) + && + current_grid_id == DataSetGetObjectID( pgpdt, this->internalgrid_indices[i]) ) { + + /* great, still valid! since we have the index position in + gpdt, just use it to get, save, the grid pointer: */ + + this->pg[i] = DataSetGetObjectPtr( pgpdt, this->internalgrid_indices[i]); + } + } + if ( UNDEF != (int)this->internalgrid_indicesb[i] ) { + + /* does index still point to a valid grid? (i.e., the one with + the same id): */ + + if ( InternalGridEnum() == DataSetGetEnum( pgpdtb, this->internalgrid_indices[i]) + && + current_grid_id == DataSetGetObjectID( pgpdtb, this->internalgrid_indices[i]) ) { + + /* great, still valid! since we have the index position in + gpdtb, just use it to get, save, the grid pointer: */ + + this->pg[i] = DataSetGetObjectPtr( pgpdtb, this->internalgrid_indicesb[i]); + } + } + + /* if we haven't yet "recovered" the InternalGrid object pointer, + go seek it out: */ + + if ( !this->pg[i] ) { + + foundgrid = DataSetSeek( (void**)(this->pg+i), /* cast because some compilers recognise as */ + /* ptr to ptr to struct otherwise */ + (this->internalgrid_indices+i) , + pgpdt, current_grid_id, 1, InternalGridEnum()); + if(!foundgrid){ + /*We did not find this grid in pgpdt, go find it in pgpdtb:*/ + foundgrid = DataSetSeek( (void**)(this->pg+i), /* cast because some compilers recognise as */ + /* ptr to ptr to struct otherwise */ + (this->internalgrid_indicesb+i) , + pgpdtb, current_grid_id, 1, InternalGridEnum()); + } + + if(!foundgrid){ + /*We could not find our grid in pgpdt and pgpdtb, this is not good!:*/ + noerr=0; + } + if (!noerr) { + if (first) { + printf(" %s error:\n", __func__); + first = 0; + } + printf(" could not find grid with id of %d for element %d of type \"%s\"\n", + current_grid_id, PenpairLoadGetID(this), PenpairLoadGetName(this)); + } + } + } + } + } + else { + + if (first) { + printf(" %s error:\n", __func__); + first=0; + } + printf(" could not find DataSets of type \"gpdt\", grid object\n"); + printf(" pointers not established for penpair element id = %d\n", + PenpairLoadGetID(this)); + noerr = 0; + } + + #endif + + /*Now look for the two elements that make up the penalty pair, if we are running numdof==2: */ + if(this->penpair.numdofs==2){ + found=0; + for ( j=0; jinternalelement_indices[i]){ + + this->elements[i] = DataSetGetObjectPtr( est, this->internalelement_indices[i]); + + } + + if ( !this->elements[i]){ + + noerr = DataSetSeek( (void**)(&this->elements[i]), //cast because some compilers recognise as ptr to ptr to struct otherwise + &this->internalelement_indices[i], est, this->penpair.eids[i], 1, TriaElementEnum()); + + } + + if (!noerr) { + printf("%s%s%i\n",__FUNCT__," error message: could not find element with id ",this->penpair.eids[i]); + goto cleanup_and_return; + } + } + } + else{ + printf(" could not find DataSet of type \"est\", element object\n"); + printf(" pointers not established for penpair load id = %d\n", PenpairLoadGetID(this)); + noerr=0;goto cleanup_and_return; + } + } + + #ifdef _DEBUG_ + printf(" load data after all lookups:\n"); + PenpairLoadEcho( this); + #endif + + xfree((void**)&arglist); + } + else { + + /* quiet return */ + ; + } + cleanup_and_return: + + return noerr; + +} + +/*-------------------------------------------------- + PenpairLoadPenaltyCreateKMatrix + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyCreateKMatrix" +int PenpairLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){ + + + PenpairLoad* this = NULL; + int noerr = 1; + + /*PenpairLoad between two grids. Two cases, numdofs=1, we just penalize the same dof on both grids. + * numdofs=2, the penalization is to enfore non penetration between the two grids. :*/ + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + if (this->penpair.numdofs==1){ + PenpairLoadPenaltyCreateKMatrixOneDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type); + } + else if (this->penpair.numdofs==2){ + PenpairLoadPenaltyCreateKMatrixTwoDof( pKe_gg,vpthis,inputs,K_flag,kmax,analysis_type); + } + else{ + printf("%s%s\n",__FUNCT__," error message: numdofs not supported yet!"); + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + PenpairLoadPenaltyCreateKMatrixOneDof + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixOneDof" +int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type){ + + + /* PenpairLoad is, as the name describes, a penalty load between two grids, for a certain degree of freedom. + * This is the same as adding a infinitely stiff spring element between two grids, for a certain degree + * of freedom. The goal is to make the value of that dof identical on both grids. + * In order to do so, we add a stiffness matrix of the type [1 -1 ; -1 1]*lambda, where lambda is a penalty + * parameter.*/ + + int i; + int noerr=1; + + PenpairLoad* this = NULL; + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int* structural_dof_list = NULL; + int dof; + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + /*Create element stiffness matrix for that penalty: */ + Ke_gg=NewElemMatrix(2); + + for (i=0;i<2;i++){ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[i]); + dof=structural_dof_list[this->penpair.dof-1]; //matlab indexing + Ke_gg->row_indices[i]=dof; + } + *(Ke_gg->terms+Ke_gg->nrows*0+0)=kmax*pow(10,this->penpair.penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*0+1)=-kmax*pow(10,this->penpair.penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*1+0)=-kmax*pow(10,this->penpair.penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*1+1)=kmax*pow(10,this->penpair.penalty_offset); + + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + return noerr; +} + +/*-------------------------------------------------- + PenpairLoadPenaltyCreateKMatrixTwoDof + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyCreateKMatrixTwoDof" +int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){ + + + int i,j; + int noerr=1; + + PenpairLoad* this = NULL; + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int* structural_dof_list = NULL; + int dof; + int grids[2]; + double* thickness_param=NULL; + double thickness; + double friction,length,penalty_offset; + double normal[2]; + TriaElement* triaelem=NULL; + int grid_number; + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + friction=this->penpair.friction; + length=this->penpair.length; + penalty_offset=this->penpair.penalty_offset; + normal[0]=this->penpair.normal[0]; + normal[1]=this->penpair.normal[1]; + if(this->active){ + /*There is contact, we need to constrain the normal velocities (zero penetration), and the + *contact slip friction. */ + + #ifdef _DEBUG_ + printf("Dealing with grid pair (%i,%i)\n",InternalGridGetID(this->pg[0]),InternalGridGetID(this->pg[1])); + #endif + + /*Recover input parameters: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + + /*Recover thickness at grid 1. We assume thickness at grid 2 is the same as thickness at grid1: */ + if(thickness_param){ + /*Go pickup the thickness in the inputs: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[0]); + dof=structural_dof_list[0]; + thickness=thickness_param[dof]; + } + else{ + + /*Figure out place of grid 1 of the pair in the element 1 of the pair: */ + triaelem=this->elements[0]; + for(i=0;i<3;i++){ + if (InternalGridGetID(this->pg[0])==InternalGridGetID(triaelem->pg[i])){ + grid_number=i; + break; + } + } + + /*Go pick up the thickness in the triaelem: */ + thickness=triaelem->tria.h[grid_number]; + } + + #ifdef _DEBUG_ + printf("Thickness at grid (%i,%i): %lg\n",InternalGridGetID(this->pg[0]),InternalGridGetID(this->pg[1]),thickness); + #endif + + /*Create element stiffness matrix: */ + Ke_gg=NewElemMatrix(2*NDOF2); //penalty on two grids + + for (i=0;i<2;i++){ + /* We use the structural dof list, from which we keep only the first two dof (x for vx, and y for vy).*/ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=dof; + } + } + + /*From Peter Wriggers book (Computational Contact Mechanics, p191): */ + //First line: + *(Ke_gg->terms+Ke_gg->nrows*0+0)+=pow(normal[0],2)*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*0+1)+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*0+2)+=-pow(normal[0],2)*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*0+3)+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + //Second line: + *(Ke_gg->terms+Ke_gg->nrows*1+0)+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*1+1)+=pow(normal[1],2)*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*1+2)+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*1+3)+=-pow(normal[1],2)*kmax*pow(10,penalty_offset); + //Third line: + *(Ke_gg->terms+Ke_gg->nrows*2+0)+=-pow(normal[0],2)*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*2+1)+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*2+2)+=pow(normal[0],2)*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*2+3)+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + //Fourth line: + *(Ke_gg->terms+Ke_gg->nrows*3+0)+=-normal[0]*normal[1]*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*3+1)+=-pow(normal[1],2)*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*3+2)+=normal[0]*normal[1]*kmax*pow(10,penalty_offset); + *(Ke_gg->terms+Ke_gg->nrows*3+3)+=pow(normal[1],2)*kmax*pow(10,penalty_offset); + + /*Now take care of the friction: of type sigma=friction*(tangent_velocity2-tangent_velocity1)*/ + + //First line: + *(Ke_gg->terms+Ke_gg->nrows*0+0)+=pow(normal[1],2)*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*0+1)+=-normal[0]*normal[1]*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*0+2)+=-pow(normal[1],2)*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*0+3)+=normal[0]*normal[1]*thickness*length*friction; + //Second line: + *(Ke_gg->terms+Ke_gg->nrows*1+0)+=-normal[0]*normal[1]*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*1+1)+=pow(normal[0],2)*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*1+2)+=normal[0]*normal[1]*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*1+3)+=-pow(normal[0],2)*thickness*length*friction; + //Third line: + *(Ke_gg->terms+Ke_gg->nrows*2+0)+=-pow(normal[1],2)*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*2+1)+=normal[0]*normal[1]*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*2+2)+=pow(normal[1],2)*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*2+3)+=-normal[0]*normal[1]*thickness*length*friction; + //Fourth line: + *(Ke_gg->terms+Ke_gg->nrows*3+0)+=normal[0]*normal[1]*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*3+1)+=-pow(normal[0],2)*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*3+2)+=-normal[0]*normal[1]*thickness*length*friction; + *(Ke_gg->terms+Ke_gg->nrows*3+3)+=pow(normal[0],2)*thickness*length*friction; + + } + else{ + /*the grids on both sides of the rift do not penetrate. PenpairLoadPenaltyCreatePVectorTwoDof will + *take care of adding point loads to simulate pressure on the rift flanks. But as far as stiffness, + there is none (0 stiffness implies decoupling of the flank rifts, which is exactly what we want): */ + *pKe_gg=NULL; + } + + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + return noerr; +} + + +/*-------------------------------------------------- + PenpairLoadPenaltyCreatePVector + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyCreatePVector" + +int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ + + PenpairLoad* this = NULL; + int noerr = 1; + + /*PenpairLoad between two grids. Two cases, numdofs=1, do nothing, + * numdofs=2, the penalization is to enfore non penetration between the two grids. :*/ + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + if (this->penpair.numdofs==1){ + *ppe_g=NULL; + } + else if (this->penpair.numdofs==2){ + PenpairLoadPenaltyCreatePVectorTwoDof( ppe_g,vpthis,inputs,kmax,analysis_type); + } + else{ + printf("%s%s\n",__FUNCT__," error message: numdofs not supported yet!"); + *ppe_g=NULL; + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + PenpairLoadPenaltyCreatePVectorTwoDof + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyCreatePVectorTwoDof" + +int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ + + int i,j; + int noerr=1; + + PenpairLoad* this = NULL; + Matpar* matpar = NULL; + TriaElement* triaelem=NULL; + + + /* output: */ + ElemVector* pe_g = NULL; + + int* structural_dof_list = NULL; + double* thickness_param=NULL; + double thickness; + double* bed_param=NULL; + double bed; + int dof; + int grids[2]; + double length; + double normal[2]; + int fill; + double rho_ice,gravity,rho_water; + int grid_number; + double pressure; + + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + triaelem=this->elements[0]; + matpar=triaelem->matpar; + + if(!this->active){ + /*Ok, this rift is opening. We should put loads on both sides of the rift flanks. Because we are dealing with contact mechanics, + * and we want to avoid zigzagging of the loads, we want lump the loads onto grids, not onto surfaces between grids.:*/ + + #ifdef _DEBUG_ + _printf_("Grids (%i,%i) are free of constraints\n",InternalGridGetID(this->pg[0]),InternalGridGetID(this->pg[1])); + #endif + + /*Ok, to compute the pressure, we are going to need material properties, thickness and bed for the two grids. We assume those properties to + * be the same across the rift.: */ + + length=this->penpair.length; + normal[0]=this->penpair.normal[0]; + normal[1]=this->penpair.normal[1]; + fill=this->penpair.fill; + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + gravity=matpar->g; + + /*get thickness: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + if(thickness_param){ + /*Go pickup the thickness in the inputs: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[0]); + dof=structural_dof_list[0]; + thickness=thickness_param[dof]; + } + else{ + /*Figure out place of grid 1 of the pair in the element 1 of the pair: */ + for(i=0;i<3;i++){ + if (InternalGridGetID(this->pg[0])==InternalGridGetID(triaelem->pg[i])){ + grid_number=i; + break; + } + } + /*Go pick up the thickness in the triaelem: */ + thickness=triaelem->tria.h[grid_number]; + } + + /*get bed: */ + bed_param=ParameterInputsRecover(inputs,"bed"); + if(bed_param){ + /*Go pickup the bed in the inputs: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[0]); + dof=structural_dof_list[0]; + bed=bed_param[dof]; + } + else{ + /*Figure out place of grid 1 of the pair in the element 1 of the pair: */ + for(i=0;i<3;i++){ + if (InternalGridGetID(this->pg[0])==InternalGridGetID(triaelem->pg[i])){ + grid_number=i; + break; + } + } + /*Go pick up the bed in the triaelem: */ + bed=triaelem->tria.h[grid_number]; + } + + + /*Ok, now compute the pressure (in norm) that is being applied to the flanks, depending on the type of fill: */ + if(fill==WATERFILL){ + if(triaelem->tria.shelf){ + /*We are on an ice shelf, hydrostatic equilibrium is used to determine the pressure for water fill: */ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(bed,(double)2)/(double)2; + } + else{ + //We are on an icesheet, we assume the water column fills the entire front: */ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2 - rho_water*gravity*pow(thickness,(double)2)/(double)2; + } + } + else if(fill==AIRFILL){ + pressure=rho_ice*gravity*pow(thickness,(double)2)/(double)2; //icefront on an ice sheet, pressure imbalance ice vs air. + } + else if(fill==ICEFILL){ //icefront finding itself against another icefront (pressure imbalance is fully compensated, ice vs ice) + pressure=0; + } + else{ + _printf_("%s%s%i%s\n",__FUNCT__," error message: fill type ",fill," not supported yet."); + noerr=0;goto cleanup_and_return; + } + + /*Ok, now we have the pressure and all the properties, create the vector, with correct partitioning of degrees of freedom: */ + pe_g=NewElemVector(4); //2 grids, 2 dofs in vx and vy + for (i=0;i<2;i++){ + structural_dof_list= InternalGridGetStructuralDoflistPtr( this->pg[i]); + for (j=0;j<2;j++){ + dof=structural_dof_list[j]; + *(pe_g->row_indices+i*2+j)=dof; + } + } + + /*Ok, add contribution to first grid, along the normal i==0: */ + for (j=0;j<2;j++){ + *(pe_g->terms+0*2+j)+=pressure*normal[j]*length; + } + + /*Add contribution to second grid, along the opposite normal: i==1 */ + for (j=0;j<2;j++){ + *(pe_g->terms+1*2+j)+= -pressure*normal[j]*length; + } + } + else{ + /*The penalty is active. No loads implied here.*/ + pe_g=NULL; + } + cleanup_and_return: + + /*Assign output pointer:*/ + *ppe_g=pe_g; + return noerr; +} +int PotentialUnstableConstraints(DataSet* lst,ParameterInputs* inputs,int analysis_type_enum); + +/*-------------------------------------------------- + PenpairLoadPotentialUnstableConstraint + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPotentialUnstableConstraint" + +int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ + + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = NULL; + + int noerr=1; + + /*parameters: */ + double* velocity=NULL; + + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + int i,j; + int free_boundaries; + double vxvy_list[2][2]; //velocities for all grids + double normal[2]; + int dof; + double vx1,vy1; + double vx2,vy2; + double penetration; + int* structural_dof_list = NULL; + + int unstable; + int active; + + /*First recover parameter inputs: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + + if(!velocity){ + printf("%s%s\n",__FUNCT__," error message: missing velocity input!"); + return 0; + } + + normal[0]=this->penpair.normal[0]; + normal[1]=this->penpair.normal[1]; + + /*Recover velocities for both grids: */ + for (i=0;i<2;i++){ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[i]); + for (j=0;jpg[i]),vxvy_list[i][0],vxvy_list[i][1]); + } + printf("Normal for grids: [%g,%g]\n",normal[0],normal[1]); + #endif + + vx1=vxvy_list[0][0]; vy1=vxvy_list[0][1]; + vx2=vxvy_list[1][0]; vy2=vxvy_list[1][1]; + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1]; + + /*Ok, we are looking for positive penetration in an active constraint: */ + if(this->active){ + if (penetration>=0){ + unstable=1; + } + else{ + unstable=0; + } + } + else{ + unstable=0; + } + + cleanup_and_return: + /*assign output pointer: */ + *punstable=unstable; + + return noerr; +} +/*-------------------------------------------------- + PenpairLoadPenaltyPreConstrain + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyPreConstrain" + +int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ + + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = NULL; + + int noerr=1; + + /*parameters: */ + double* velocity=NULL; + + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + int i,j; + int free_boundaries; + double vxvy_list[2][2]; //velocities for all grids + double normal[2]; + int dof; + double vx1,vy1; + double vx2,vy2; + double penetration; + int* structural_dof_list = NULL; + + int unstable; + int active; + + /*First recover parameter inputs: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + + if(!velocity){ + printf("%s%s\n",__FUNCT__," error message: missing velocity input!"); + return 0; + } + + normal[0]=this->penpair.normal[0]; + normal[1]=this->penpair.normal[1]; + + /*Recover velocities for both grids: */ + for (i=0;i<2;i++){ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[i]); + for (j=0;jpg[i]),vxvy_list[i][0],vxvy_list[i][1]); + } + printf("Normal for grids: [%g,%g]\n",normal[0],normal[1]); + #endif + + vx1=vxvy_list[0][0]; vy1=vxvy_list[0][1]; + vx2=vxvy_list[1][0]; vy2=vxvy_list[1][1]; + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1]; + + + /*Ok, we are preconstraining here. Ie, anything that penetrates is constrained until stability of the entire set + * of constraints is reached.: */ + if(penetration<0){ + if (!this->active){ + /*This is the first time penetration happens: */ + this->active=1; + unstable=1; + } + else{ + /*This constraint was already active: */ + this->active=1; + unstable=0; + } + } + else{ + /*No penetration happening. : */ + if (!this->active){ + /*This penalty was not active, and no penetration happening. Do nonthing: */ + this->active=0; + unstable=0; + } + else{ + /*Ok, this penalty wants to get released. But not now, this is preconstraint, not constraint: */ + this->active=1; + unstable=0; + } + } + + cleanup_and_return: + /*assign output pointer: */ + *punstable=unstable; + + return noerr; +} + + +/*-------------------------------------------------- + PenpairLoadPenaltyConstrain + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenaltyConstrain" + +int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ + + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = NULL; + + int noerr=1; + + /*parameters: */ + double* velocity=NULL; + + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + int i,j; + int free_boundaries; + double vxvy_list[2][2]; //velocities for all grids + double normal[2]; + int dof; + double vx1,vy1; + double vx2,vy2; + double penetration; + int* structural_dof_list = NULL; + + int unstable; + int active; + double* max_penetrationparam; + double max_penetration; + + /*First recover parameter inputs: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + max_penetrationparam=ParameterInputsRecover(inputs,"max_penetration"); + + if(!velocity){ + printf("%s%s\n",__FUNCT__," error message: missing velocity input!"); + return 0; + } + if(!max_penetrationparam){ + printf("%s%s\n",__FUNCT__," error message: missing max_penetration input!"); + return 0; + } + max_penetration=*max_penetrationparam; + + normal[0]=this->penpair.normal[0]; + normal[1]=this->penpair.normal[1]; + + /*Recover velocities for both grids: */ + for (i=0;i<2;i++){ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[i]); + for (j=0;jpg[i]),vxvy_list[i][0],vxvy_list[i][1]); + } + printf("Normal for grids: [%g,%g]\n",normal[0],normal[1]); + #endif + + vx1=vxvy_list[0][0]; vy1=vxvy_list[0][1]; + vx2=vxvy_list[1][0]; vy2=vxvy_list[1][1]; + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1]; + + /*Activate or deactivate penalties: */ + if(penetration<0){ + /*There is penetration, we need to active the penalty so this penetration will be NULL: */ + active=1; + } + else{ + /*Ok, there is no penetration for these two grids of the rift. We want to deactivate the penalty. If we do + * it all at once, then we will zigzag forever. Only deactivate once at a time. Deactivate the one that most + * wants to, ie the one with the max penetration: */ + if (penetration==max_penetration){ + active=0; + } + else{ + /*Ok, we are dealing with the rest of the penalties that want to be freed. But may be in this lot there + * are penalties that were already free? Keep them as such. For the rest, do not release them yet: */ + if (this->active==0){ + active=0; + } + else{ + active=1; + } + } + } + + /*Here, we try to avoid zigzaging. When a penalty activates and deactivates for more than penalty_lock times, + * we lock it: */ + if(this->counter>this->penpair.penalty_lock){ + /*This penalty has zig zagged too long, fix it to smooth results: */ + active=1; this->active=active; + printf(" locking penalty %i\n",PenpairLoadGetID(this)); + } + + //Figure out stability of this penalty + if(this->active==active){ + unstable=0; + } + else{ + unstable=1; + this->counter++; + } + + //Set penalty flag + this->active=active; + + + cleanup_and_return: + /*assign output pointer: */ + *punstable=unstable; + + return noerr; +} + +/*-------------------------------------------------- + PenpairLoadPenetration + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadPenetration" + +int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type){ + + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = NULL; + + int noerr=1; + + /*parameters: */ + double* velocity=NULL; + + + /*Some pointer intialization: */ + this = (PenpairLoad*)vpthis; + + int i,j; + int free_boundaries; + double vxvy_list[2][2]; //velocities for all grids + double normal[2]; + int dof; + double vx1,vy1; + double vx2,vy2; + double penetration; + int* structural_dof_list = NULL; + + int unstable; + int active; + + /*First recover parameter inputs: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + if(velocity){ + normal[0]=this->penpair.normal[0]; + normal[1]=this->penpair.normal[1]; + + /*Recover velocities for both grids: */ + for (i=0;i<2;i++){ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg[i]); + for (j=0;jpg[i]),vxvy_list[i][0],vxvy_list[i][1]); + } + printf("Normal for grids: [%g,%g]\n",normal[0],normal[1]); + #endif + + vx1=vxvy_list[0][0]; vy1=vxvy_list[0][1]; + vx2=vxvy_list[1][0]; vy2=vxvy_list[1][1]; + + /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ + penetration=(vx2-vx1)*normal[0]+(vy2-vy1)*normal[1]; + + /*Now, we return penetration only if we are active!: */ + if(this->active==0)penetration=0; + } + else{ + printf("%s%s\n",__FUNCT__," error message: missing velocity input!"); + return 0; + } + + + cleanup_and_return: + /*assign output pointer: */ + *ppenetration=penetration; + + return noerr; +} +/*-------------------------------------------------- + PenpairLoadGetName + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PenpairLoadGetName" + +char* PenpairLoadGetName( void* vpthis){ + + /* vpthis for polymorphic function compatibility */ + PenpairLoad* this = (PenpairLoad*)vpthis; + + return PenpairGetName(&this->penpair); + +} + +#undef NDOF2