Index: /issm/trunk/src/c/objects/cielo/IcefrontLoad.c =================================================================== --- /issm/trunk/src/c/objects/cielo/IcefrontLoad.c (revision 216) +++ /issm/trunk/src/c/objects/cielo/IcefrontLoad.c (revision 216) @@ -0,0 +1,1556 @@ + +/* + IcefrontLoad.c: +*/ + +/* IcefrontLoad: this object deals with boundary conditions at the icefront, + * for each type of runs: dynamic water pressure for diagnostic runs, free or + * set radiation for prognostic runs. */ + +/*For debugging purposes: */ +#define ELID 41 //element for which to print debug statements +#define RANK 2 //rank of cpu for which to print debug statements. +//#define _DEBUGELEMENTS_ +//#define _DEBUGGAUSS_ + + +#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 "./IcefrontLoad.h" + +/*Objects: */ +#include "../objects.h" + +#include "../ElementObjects/GetElementGridData.h" +#include "../ElementObjects/GaussPoints.h" + +#define NDOF1 1 +#define NDOF2 2 + + +/*-------------------------------------------------- + NewIcefrontLoad + --------------------------------------------------*/ + +IcefrontLoad* NewIcefrontLoad( void* vpicefront) +{ + + Icefront* icefront = (Icefront*)vpicefront; + IcefrontLoad* this = NULL; + + int i; + + this=(IcefrontLoad*)xcalloc(1,sizeof(IcefrontLoad)); /* guaranteed NULL, 0 initialisation */ + + /* include the "base class" in the "derived" object: */ + memcpy( &this->icefront, icefront, sizeof(Icefront)); + + strcpy(this->name,"IcefrontLoad"); + for(i=0;iinternalgrid_indices[i]=UNDEF; + this->internalgrid_indicesb[i]=UNDEF; + this->pg[i]=NULL; + this->element_index=UNDEF; + this->element_enum=UNDEF; + this->element=NULL; + this->matpar_index=UNDEF; + this->matpar_enum=UNDEF; + this->matpar=NULL; +} + + return this; +} + +/*-------------------------------------------------- + DeleteIcefrontLoad + --------------------------------------------------*/ + +void DeleteIcefrontLoad(void* *vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + IcefrontLoad* *pthis = (IcefrontLoad**)vpthis; + IcefrontLoad* this = *pthis; + + xfree(vpthis); +} + + +/*-------------------------------------------------- + IcefrontLoadEcho + --------------------------------------------------*/ + +void IcefrontLoadEcho(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + IcefrontLoad* this = (IcefrontLoad*)vpthis; + + int i,j; + + printf("\nIcefrontLoad echo:\n"); + printf("-------------------\n"); + printf(" base Icefront object:\n"); + IcefrontEcho( &this->icefront); + printf("\n"); + + printf(" name: %s\n" ,this->name); + printf(" load point to following grids (pointer-based ops):\n"); + printf(" local grid number corresponding grid id\n"); + printf(" ----------------- ---------------------\n"); + for ( i=0; ipg[i]) + printf(" %8d %8d\n", i, InternalGridGetID( this->pg[i])); + } + if(this->element)printf(" load element: %p enum %i\n",this->element,this->element_enum); + if(this->matpar)printf(" load matpar: %p enum %i\n",this->matpar,this->matpar_enum); + + return; +} + +/*-------------------------------------------------- + IcefrontLoadSizeOfInDoubles + --------------------------------------------------*/ + +int IcefrontLoadSizeOfInDoubles(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + IcefrontLoad* this = (IcefrontLoad*)vpthis; + int struct_doubles=0; + + struct_doubles = + (int)ceil((float)sizeof(*this)/(float)sizeof(double)); + + #ifdef _DEBUG_ + printf("\n IcefrontLoadSizeOfInDoubles diagnostics:\n"); + printf(" returning %d doubles\n", struct_doubles); + #endif + + return struct_doubles; +} + + +/*-------------------------------------------------- + IcefrontLoadShrinkWrap + --------------------------------------------------*/ + +void IcefrontLoadShrinkWrap(void* to, void* vpthis) +{ + /* no "extended" data means a simple copy will do: */ + + memcpy( to, vpthis, sizeof(double)*IcefrontLoadSizeOfInDoubles(vpthis)); +} + + +/*-------------------------------------------------- + IcefrontLoadUnwrap + --------------------------------------------------*/ + +int IcefrontLoadUnwrap(void* vpthis) +{ + /* no extended data means no internal pointers to reestablish */ + + return 1; +} + +/*-------------------------------------------------- + IcefrontLoadMarshall + --------------------------------------------------*/ + +void IcefrontLoadMarshall( char* *pbuffer, void* vpthis, int* size) +{ + /* vpthis for polymorphic function compatibility */ + IcefrontLoad* this = (IcefrontLoad*)vpthis; + + char* buffer=*pbuffer; + + int i, j, size_icefront=0; + int size_string1; + int ig; + + int num_int=0, num_double=0; + + /* first marshall the embedded Icefront object: ... */ + IcefrontMarshall( &buffer, &this->icefront, &size_icefront); + + /* ... then marshall the IcefrontLoad member data: */ + + size_string1=strlen(this->name)+1; + cellmar(&size_string1,&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->name,&buffer,STRING_FLAG); + for(i=0;iinternalgrid_indices[i],&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->internalgrid_indicesb[i],&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->pg[i],&buffer,INTEGER_FLAG);num_int++; + } + cellmar(&this->element_index,&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->element_enum,&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->element,&buffer,INTEGER_FLAG);num_int++; + + cellmar(&this->matpar_index,&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->matpar_enum,&buffer,INTEGER_FLAG);num_int++; + cellmar(&this->matpar,&buffer,INTEGER_FLAG);num_int++; + + /* return values: */ + *size=size_string1*sizeof(char)+ num_int*sizeof(int)+num_double*sizeof(double); + *pbuffer=buffer; + + return; +} + +/*-------------------------------------------------- + IcefrontLoadDemarshall + --------------------------------------------------*/ + +void* IcefrontLoadDemarshall( DataSet* ds, char* *pbuffer, int machine_flag) +{ + int i, j, ig; + int size_string1; + Icefront* picefront=NULL; /* for embedded object */ + + char* buffer=*pbuffer; + IcefrontLoad* this=NULL; + + /* first demarshall the embedded Icefront object: ... */ + + picefront = IcefrontDemarshall( NULL, &buffer, machine_flag); + /* create the derived object: */ + this = NewIcefrontLoad( picefront); + /* and free icefront now that it's embedded in the derived object: */ + DeleteIcefront( (void**)&picefront); + + /* ... then demarshall the IcefrontLoad member data: */ + + + celldem(&size_string1,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->name,&buffer,STRING_FLAG,machine_flag,size_string1); + for(i=0;iinternalgrid_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); + this->pg[i]=NULL; + } + celldem(&this->element_index,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->element_enum,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->element,&buffer,INTEGER_FLAG,machine_flag); + + celldem(&this->matpar_index,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->matpar_enum,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->matpar,&buffer,INTEGER_FLAG,machine_flag); + + if (ds) IcefrontLoadAddToDataSet( ds, this); + *pbuffer=buffer; + + return (void*)this; +} + +/*-------------------------------------------------- + IcefrontLoadGetID + --------------------------------------------------*/ + +int IcefrontLoadGetID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + IcefrontLoad* this = (IcefrontLoad*)vpthis; + + return IcefrontGetID(&this->icefront); +} + +/*-------------------------------------------------- + IcefrontLoadGetElementID + --------------------------------------------------*/ + +int IcefrontLoadGetElementID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + IcefrontLoad* this = (IcefrontLoad*)vpthis; + + return IcefrontGetElementID(&this->icefront); +} + + +/*-------------------------------------------------- + IcefrontLoadSetVirtualFunctions + --------------------------------------------------*/ + +int IcefrontLoadSetVirtualFunctions( void* vp) +{ + + /* VirtualFunctions implemented as a "friend" class: */ + + /* void* for polymorphic function compatibility */ + VirtualFunctions* vf = (VirtualFunctions*)vp; + + if ( VirtualFunctionsInit(vf) ) { + + /* set IcefrontLoad-specific virtual function pointers: */ + + /* general: */ + + vf->GetName = &IcefrontLoadGetName; + + /* configuration: */ + vf->Configure = &IcefrontLoadConfigure; + + /* loads generation: */ + vf->CreatePVector = &IcefrontLoadCreatePVector; + + return 1; + } + else + return 0; + +} + +/*-------------------------------------------------- + IcefrontLoadGetName + --------------------------------------------------*/ + +char* IcefrontLoadGetName( void* vpthis) +{ + IcefrontLoad* this = (IcefrontLoad*)vpthis; + + #ifdef _DEBUG_ + printf("virtual function: IcefrontLoadGetName\n"); + #endif + + return this->name; +} + + +/*-------------------------------------------------- + IcefrontLoadAddToDataSet + --------------------------------------------------*/ + +int IcefrontLoadAddToDataSet( DataSet* dataset, IcefrontLoad* 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, IcefrontLoadEnum(), (void*)this, + &DeleteIcefrontLoad, &IcefrontLoadEcho, &IcefrontLoadSizeOfInDoubles, + &IcefrontLoadShrinkWrap, &IcefrontLoadUnwrap, + &IcefrontLoadMarshall, + &IcefrontLoadGetID, &IcefrontLoadSetVirtualFunctions); + + return 1; +} + + +/*-------------------------------------------------- + IcefrontLoadConfigure + --------------------------------------------------*/ + +int IcefrontLoadConfigure( void* vpthis, int num_datasets, ...) { + + /* set up data connections to all necessary, related objects: */ + + IcefrontLoad* this = (IcefrontLoad*)vpthis; + + int i,j, found, foundgrid, noerr=1, first=1; + + /* expected DataSet input: */ + + DataSet* pgpdt=NULL; + DataSet* pgpdtb=NULL; + DataSet* pest=NULL; + DataSet* pmpt=NULL; + + /* variable argument list handling: */ + + va_list args; + DataSet* *arglist=NULL; + + #ifdef _DEBUG_ + printf("virtual function: IcefrontLoadConfigure\n"); + #endif + + + if (this) { + + if (!num_datasets) { + + printf(" %s error:\n", __func__); + 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;ipg[i]=NULL; + } + + /* grid connections based on DataSet of type "gpdt": */ + + #ifndef _PARALLEL_ + found=0; + for ( i=0; iicefront.g; + + for ( i=0; iicefront); 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) { + if (first) { + printf(" %s error:\n", __func__); + first = 0; + } + printf(" could not find grid with id of %d for load %d of type \"%s\"\n", + current_grid_id, IcefrontLoadGetID(this), IcefrontLoadGetName(this)); + } + } + } + } + } + else { + + if (first) { + printf(" %s error:\n", __func__); + first=0; + } + printf(" could not find DataSet of type \"gpdt\", grid object\n"); + printf(" pointers not established for penpair load id = %d\n", + IcefrontLoadGetID(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; iicefront.g; + + for ( i=0; iicefront); 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_indices[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, IcefrontLoadGetElementID(this), IcefrontLoadGetName(this)); + if(current_grid_id==306)IcefrontLoadEcho(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", + IcefrontLoadGetID(this)); + noerr = 0; + } + + #endif + + + /*Now look for the element pointer corresponding to this load: */ + found=0; + for ( j=0; jelement_index) { + this->element = DataSetGetObjectPtr( pest, this->element_index); + this->element_enum = DataSetGetEnum( pest, this->element_index); + } + + if ( !this->element) { + + noerr *= DataSetSeek( (void**)(&this->element), //cast because some compilers recognise as ptr to ptr to struct otherwise + &this->element_index , pest, IcefrontLoadGetElementID(this), 0); + this->element_enum = DataSetGetEnum( pest, this->element_index); + } + if (!noerr) { + printf(" Icefront load %i could not find element with id %d \n", IcefrontLoadGetID(this),IcefrontLoadGetElementID(this)); + goto cleanup_and_return; + } + + /* matpar property connections based on DataSet of type "mpt": */ + found=0; + for ( i=0; iicefront); + + if ( this->matpar_index!= UNDEF){ + + + if ( this->matpar_enum == DataSetGetEnum( pmpt, this->matpar_index) + && + current_matpar_id == DataSetGetObjectID( pmpt, this->matpar_index) ) { + + /* index is still valid; just get matpar entry pointer directly: */ + + this->matpar = DataSetGetObjectPtr( pmpt, this->matpar_index); + } + } + + /* if we haven't yet "recovered" the matpar object pointer, go seek it out: */ + + if ( NULL == this->matpar) { + + noerr *= DataSetSeek( &this->matpar, &this->matpar_index, + pmpt, current_matpar_id, 1, MatparEnum() ); /* 0 to generalise later to Pcomp, etc. ... */ + + if (!noerr) { + if (first) { + _printf_(" %s error:\n", __func__); + first = 0; + } + _printf_(" could not find matpar entry with id of %d for element %d of type \"%s\"\n", + current_matpar_id, IcefrontGetID(&this->icefront), IcefrontGetName(&this->icefront)); + } + else { + + /* though the answer for now is just MatparEnum(), for generality go ahead and + ask the question, "what type of matpar has been referenced?" */ + + this->matpar_enum = DataSetGetEnum( *(arglist+i), this->matpar_index); + } + } + + } + else { + + if (first) { + _printf_(" %s error:\n", __func__); + first=0; + } + _printf_(" could not find DataSet of type \"mpt\"; matice and matpar relationships\n"); + _printf_(" not established for penta element id = %d\n", IcefrontGetElementID(this)); + noerr = 0; + } + + #ifdef _DEBUG_ + printf(" load data after all lookups:\n"); + IcefrontLoadEcho( this); + #endif + + xfree((void**)&arglist); + + } + else { + + /* quiet return */ + ; + } + cleanup_and_return: + + return noerr; + +} + +/*-------------------------------------------------- + IcefrontLoadCreatePVector + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "IcefrontLoadCreatePVector" + +int IcefrontLoadCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){ + int noerr=1; + + IcefrontLoad* load=(IcefrontLoad*)vpthis; + //printf("Eid for icefront load: %i rank %i\n",load->icefront.eid,my_rank); + + /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */ + if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ + noerr=IcefrontLoadCreatePVectorDiagnosticHoriz( ppe_g, vpthis, inputs); + } + else{ + printf("%s%s%i%s\n",__FUNCT__," error message: analysis ",analysis_type," not supported yet"); + noerr=0; + } + return noerr; + +} + +/*-------------------------------------------------- + IcefrontLoadCreatePVectorDiagnosticHoriz + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "IcefrontLoadCreatePVectorDiagnosticHoriz" + +int IcefrontLoadCreatePVectorDiagnosticHoriz( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + + int noerr=1; + IcefrontLoad* this=NULL; + + /*First recover pointer to icefront: */ + this = (IcefrontLoad*)vpthis; + + /*Branck on the type of icefront: */ + if strcasecmp_eq(this->icefront.type,"segment"){ + noerr=IcefrontLoadCreatePVectorDiagnosticHorizSegment(ppe_g,vpthis, inputs); + } + else{ + noerr=IcefrontLoadCreatePVectorDiagnosticHorizQuad(ppe_g,vpthis, inputs); + } + return noerr; +} + +/*-------------------------------------------------- + IcefrontLoadCreatePVectorDiagnosticHorizSegment + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "IcefrontLoadCreatePVectorDiagnosticHorizSegment" + +int IcefrontLoadCreatePVectorDiagnosticHorizSegment( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + int noerr=1; + int i,j; + + int numgrids=2; + int numdof=numgrids*NDOF2; + + int* structural_dof_list = NULL; + int dof; + + /* output: */ + ElemVector* pe_g = NULL; + + IcefrontLoad* this=NULL; + Matpar* matpar=NULL; + TriaElement* triaelement=NULL; + PentaElement* pentaelement=NULL; + + /* material parameters: */ + double rho_ice; + double rho_water; + double gravity; + + double x1,y1,x2,y2; + + /* input parameters: */ + double* thickness_param=NULL; + double thickness_list_element[3]; + double thickness_list[2]; + double* bed_param=NULL; + double bed_list_element[3]; + double bed_list[2]; + int grid1,grid2; + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + + double normal[2]; + double length; + + int fill; + + /*First recover pointer to icefront: */ + this = (IcefrontLoad*)vpthis; + + /*Recover material and fill parameters, according to which element the icefront belongs to: */ + if (this->icefront.element_type==TriaEnum()){ + triaelement=this->element; + matpar=triaelement->matpar; + fill=triaelement->tria.shelf; + } + else{ + pentaelement=this->element; + matpar=pentaelement->matpar; + fill=pentaelement->penta.shelf; + } + + //check that the element is onbed (collapsed formulation) otherwise:pe=0 + if(this->icefront.element_type==PentaEnum()){ + if (!pentaelement->penta.onbed){ + *ppe_g=NULL; + return noerr; + } + } + + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + gravity=matpar->g; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /*recover extra inputs: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + bed_param=ParameterInputsRecover(inputs,"bed"); + + /* Get all element grid data */ + GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + + /*Identify which grids are comprised in the segment. Watch out, different logic if we are in a penta, or a tria: */ + if (this->icefront.element_type==TriaEnum()){ + /*go through all grids and identify 1st and 2nd grid: */ + for(i=0;i<3;i++){ + if (triaelement->pg[i]==this->pg[0])grid1=i; + if (triaelement->pg[i]==this->pg[1])grid2=i; + } + } + else{ + /*go through first 3 grids, which are at the base of the penta, and identify 1st and 2nd grid: */ + for(i=0;i<3;i++){ + if (pentaelement->pg[i]==this->pg[0])grid1=i; + if (pentaelement->pg[i]==this->pg[1])grid2=i; + } + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=dof; + } + } + /*Now build thickness_list_element and bed_list_element: */ + for (i=0;i<3;i++){ + + if (this->icefront.element_type==TriaEnum()) structural_dof_list= InternalGridGetStructuralDoflistPtr( triaelement->pg[i]); + else structural_dof_list= InternalGridGetStructuralDoflistPtr( pentaelement->pg[i]); + + dof=structural_dof_list[0]; + if(thickness_param){ + thickness_list_element[i]=*(thickness_param+dof); + } + else{ + if (this->icefront.element_type==TriaEnum()) thickness_list_element[i]=triaelement->tria.h[i]; + else thickness_list_element[i]=pentaelement->penta.h[i]; + } + if(bed_param){ + bed_list_element[i]=*(bed_param+dof); + } + else{ + if (this->icefront.element_type==TriaEnum()) bed_list_element[i]= triaelement->tria.b[i]; + else bed_list_element[i]= pentaelement->penta.b[i]; + } + } + /*Build thickness_list and bed_list: */ + thickness_list[0]=thickness_list_element[grid1]; + thickness_list[1]=thickness_list_element[grid2]; + bed_list[0]=bed_list_element[grid1]; + bed_list[1]=bed_list_element[grid2]; + + //Recover grid coordinates + x1=xyz_list[0][0]; + y1=xyz_list[0][1]; + x2=xyz_list[1][0]; + y2=xyz_list[1][1]; + + //Compute length and normal of segment + normal[0]=cos(atan2(x1-x2,y2-y1)); + normal[1]=sin(atan2(x1-x2,y2-y1)); + length=sqrt(pow(x2-x1,2)+pow(y2-y1,2)); + + //Compute load contribution for this segment: + noerr=SegmentPressureLoad(pe_g,rho_water,rho_ice,gravity,thickness_list,bed_list,normal,length,fill); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && this->icefront.eid==ELID){ + printf("\nicefront load\n"); + printf("grids %i %i\n",grid1,grid2); + printf("rho_water %g\n",rho_water); + printf("rho_ice %g\n",rho_ice); + printf("gravity %g\n",gravity); + printf("thickness_list (%g,%g)\n",thickness_list[0],thickness_list[1]); + printf("bed_list (%g,%g)\n",bed_list[0],bed_list[1]); + printf("normal (%g,%g)\n",normal[0],normal[1]); + printf("length %g\n",length); + printf("fill %i\n",fill); + printf("pe_g->terms\n"); + for(i=0;iterms+i)); + } + } + #endif + + cleanup_and_return: + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build icefront load vector."); + return 0; + } + else{ + /*Assign output pointer: */ + *ppe_g=pe_g; + } + return noerr; +} + +/*-------------------------------------------------- + SegmentPressureLoad + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "SegmentPressureLoad" + +int SegmentPressureLoad(ElemVector* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill){ + int noerr=1; + + double nx,ny; + double h1,h2,b1,b2; + int num_gauss; + double* segment_gauss_coord=NULL; + double* gauss_weights=NULL; + int ig; + double Jdet; + double thickness,bed; + double ice_pressure,water_pressure,air_pressure; + double surface_under_water,base_under_water; + double pressure; + + double pe_g_gaussian[4]; + + nx=normal[0]; + ny=normal[1]; + + //Get gaussian points and weights. order 2 since we have a product of 2 nodal functions + num_gauss=3; + GaussSegment(&segment_gauss_coord, &gauss_weights, num_gauss); + + //recover thickness and bed at two grids + h1=thickness_list[0]; + h2=thickness_list[1]; + b1=bed_list[0]; + b2=bed_list[1]; + + //compute Jacobian of segment + Jdet=1./2.*length; + + for(ig=0;igterms+2*0+0)+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*nx; + *(pe_g->terms+2*0+1)+= pressure*Jdet*gauss_weights[ig]*(1+segment_gauss_coord[ig])/2*ny; + *(pe_g->terms+2*1+0)+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*nx; + *(pe_g->terms+2*1+1)+= pressure*Jdet*gauss_weights[ig]*(1-segment_gauss_coord[ig])/2*ny; + + } //for(ig=0;igicefront.element_type==TriaEnum()){ + _printf_("%s%s\n",__FUNCT__," error message: quad icefront is associated to a TriaElem!"); + return 0; + } + + /*Recover material and fill parameters: */ + pentaelement=this->element; + matpar=pentaelement->matpar; + fill=pentaelement->penta.shelf; + + /*Recover material parameters: */ + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + gravity=matpar->g; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numgridsload*NDOF2); + + /*recover extra inputs: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + bed_param=ParameterInputsRecover(inputs,"bed"); + + /* Get all element grid data */ + GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, pentaelement->pg, numgridselement, 0 ); + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=dof; + } + } + + for (i=0;ipg[i]); + dof=structural_dof_list[0]; + if(thickness_param){ + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=pentaelement->penta.h[i]; + } + if(bed_param){ + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]= pentaelement->penta.b[i]; + } + } + + //Identify which grids are comprised in the quad: + grid1=-1; grid2=-1; grid3=-1; grid4=-1; + for(i=0;ipg[i]==this->pg[0])grid1=i; + if (pentaelement->pg[i]==this->pg[1])grid2=i; + if (pentaelement->pg[i]==this->pg[2])grid3=i; + if (pentaelement->pg[i]==this->pg[3])grid4=i; + } + if((grid1==-1) || (grid2==-1)|| (grid3==-1)||(grid4==-1)){ + printf("%s%s\n",__FUNCT__," error message: could not find element grids corresponding to quad icefront!"); + return 0; + } + + /*Build new xyz, bed and thickness lists for grids 1 to 4: */ + xyz_list_quad[0][0]=xyz_list[grid1][0]; + xyz_list_quad[0][1]=xyz_list[grid1][1]; + xyz_list_quad[0][2]=xyz_list[grid1][2]; + xyz_list_quad[1][0]=xyz_list[grid2][0]; + xyz_list_quad[1][1]=xyz_list[grid2][1]; + xyz_list_quad[1][2]=xyz_list[grid2][2]; + xyz_list_quad[2][0]=xyz_list[grid3][0]; + xyz_list_quad[2][1]=xyz_list[grid3][1]; + xyz_list_quad[2][2]=xyz_list[grid3][2]; + xyz_list_quad[3][0]=xyz_list[grid4][0]; + xyz_list_quad[3][1]=xyz_list[grid4][1]; + xyz_list_quad[3][2]=xyz_list[grid4][2]; + + thickness_list_quad[0]=thickness_list[grid1]; + thickness_list_quad[1]=thickness_list[grid2]; + thickness_list_quad[2]=thickness_list[grid3]; + thickness_list_quad[3]=thickness_list[grid4]; + + bed_list_quad[0]=bed_list[grid1]; + bed_list_quad[1]=bed_list[grid2]; + bed_list_quad[2]=bed_list[grid3]; + bed_list_quad[3]=bed_list[grid4]; + + //Create a new grid in the midle of the quad and add it at the end of the list + xyz_list_quad[4][0] = (xyz_list_quad[0][0]+xyz_list_quad[1][0]+xyz_list_quad[2][0]+xyz_list_quad[3][0])/4.0; + xyz_list_quad[4][1] = (xyz_list_quad[0][1]+xyz_list_quad[1][1]+xyz_list_quad[2][1]+xyz_list_quad[3][1])/4.0; + xyz_list_quad[4][2] = (xyz_list_quad[0][2]+xyz_list_quad[1][2]+xyz_list_quad[2][2]+xyz_list_quad[3][2])/4.0; + + //Compute four normals (the quad is divided into four triangles) + v15[0]=xyz_list_quad[0][0]-xyz_list_quad[4][0]; + v15[1]=xyz_list_quad[0][1]-xyz_list_quad[4][1]; + v15[2]=xyz_list_quad[0][2]-xyz_list_quad[4][2]; + + v25[0]=xyz_list_quad[1][0]-xyz_list_quad[4][0]; + v25[1]=xyz_list_quad[1][1]-xyz_list_quad[4][1]; + v25[2]=xyz_list_quad[1][2]-xyz_list_quad[4][2]; + + v35[0]=xyz_list_quad[2][0]-xyz_list_quad[4][0]; + v35[1]=xyz_list_quad[2][1]-xyz_list_quad[4][1]; + v35[2]=xyz_list_quad[2][2]-xyz_list_quad[4][2]; + + v45[0]=xyz_list_quad[3][0]-xyz_list_quad[4][0]; + v45[1]=xyz_list_quad[3][1]-xyz_list_quad[4][1]; + v45[2]=xyz_list_quad[3][2]-xyz_list_quad[4][2]; + + cross(normal1,v15,v25); + cross(normal2,v25,v35); + cross(normal3,v35,v45); + cross(normal4,v45,v15); + + normal_norm=norm(normal1);normal1[0]=normal1[0]/normal_norm;normal1[1]=normal1[1]/normal_norm;normal1[2]=normal1[2]/normal_norm; + normal_norm=norm(normal2);normal2[0]=normal2[0]/normal_norm;normal2[1]=normal2[1]/normal_norm;normal2[2]=normal2[2]/normal_norm; + normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm; + normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm; + + + //Compute load contribution for this quad: + noerr=QuadPressureLoad(pe_g,rho_water,rho_ice,gravity,thickness_list_quad,bed_list_quad,normal1,normal2,normal3,normal4,&xyz_list_quad[0][0],fill,this->icefront.eid); + + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && this->icefront.eid==ELID){ + printf("\nicefront load\n"); + printf("grids %i %i %i %i\n",grid1,grid2,grid3,grid4); + printf("rho_water %g\n",rho_water); + printf("rho_ice %g\n",rho_ice); + printf("gravity %g\n",gravity); + printf("thickness_list (%g,%g,%g,%g)\n",thickness_list[0],thickness_list[1],thickness_list[2],thickness_list[3]); + printf("bed_list (%g,%g,%g,%g)\n",bed_list[0],bed_list[1],bed_list[2],bed_list[3]); + printf("normal1 (%g,%g,%g)\n",normal1[0],normal1[1],normal1[2]); + printf("normal2 (%g,%g,%g)\n",normal2[0],normal2[1],normal2[2]); + printf("normal3 (%g,%g,%g)\n",normal3[0],normal3[1],normal3[2]); + printf("normal4 (%g,%g,%g)\n",normal4[0],normal4[1],normal4[2]); + printf("fill %i\n",fill); + printf("pe_g->terms\n"); + for(i=0;iterms+i)); + } + } + #endif + + cleanup_and_return: + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build icefront load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + /*Assign output pointer: */ + *ppe_g=pe_g; + } + return noerr; +} + + + +/*-------------------------------------------------- + QuadPressureLoad + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "QuadPressureLoad" + +int QuadPressureLoad(ElemVector* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, + double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill,int eid){ + + + //The quad is divided into four triangles tria1 tria2 tria3 and tria4 as follows + // + // grid 4 +-----------------+ grid 3 + // |\2 1 /| + // |1\ tria3 /2| + // | \ / | + // | \ / | + // | \ / | + // | \ / | + // | \ 3 / | + // |tria4 \ / 3 | + // | 3 \grid5 | + // | / \ | + // | / 3 \ tria2| + // | / \ | + // | / \ | + // | / \ | + // | / tria1 \ | + // |2/1 2\1| + // grid1 +-----------------+ grid 2 + // + // + + int noerr=1; + int i,j; + + double nx[4]; + double ny[4]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_coord[3]; + + double surface_list[5]; + double x[5]; + double y[5]; + double z[5]; + double l12,l14,l15,l23,l25,l34,l35,l45; + double cos_theta_triangle1,cos_theta_triangle2,cos_theta_triangle3,cos_theta_triangle4; + + double xyz_tria[12][2]; + double l1l2l3_tria[3]; + double r_tria,s_tria; + double r_quad[4]; + double s_quad[4]; + double l1l4_tria[4][4]; + + double J[4]; + double z_g[4]; + double ice_pressure_tria[4]; + double air_pressure_tria; + double water_level_above_g_tria; + double water_pressure_tria; + double pressure_tria[4]; + + + + //Build the four normal vectors + nx[0]=normal1[0]; nx[1]=normal2[0]; nx[2]=normal3[0]; nx[3]=normal4[0]; + ny[0]=normal1[1]; ny[1]=normal2[1]; ny[2]=normal3[1]; ny[3]=normal4[1]; + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + //Recover the surface of the four nodes + for(i=0;i<4;i++){ + surface_list[i]=thickness_list[i]+bed_list[i]; + } + //Add surface sor the fifth point (average of the surfaces) + surface_list[4]=(surface_list[0]+surface_list[1]+surface_list[2]+surface_list[3])/4.0; + + //Recover grid coordinates + for(i=0;i<5;i++){ + x[i]=*(xyz_list+3*i+0); + y[i]=*(xyz_list+3*i+1); + z[i]=*(xyz_list+3*i+2); + } + + //Build triangles in a 2D plan before using reference elements + l12=pow((x[1]-x[0]),2)+pow((y[1]-y[0]),2)+pow((z[1]-z[0]),2); + l14=pow((x[3]-x[0]),2)+pow((y[3]-y[0]),2)+pow((z[3]-z[0]),2); + l15=pow((x[4]-x[0]),2)+pow((y[4]-y[0]),2)+pow((z[4]-z[0]),2); + l23=pow((x[2]-x[1]),2)+pow((y[2]-y[1]),2)+pow((z[2]-z[1]),2); + l25=pow((x[4]-x[1]),2)+pow((y[4]-y[1]),2)+pow((z[4]-z[1]),2); + l34=pow((x[3]-x[2]),2)+pow((y[3]-y[2]),2)+pow((z[3]-z[2]),2); + l35=pow((x[4]-x[2]),2)+pow((y[4]-y[2]),2)+pow((z[4]-z[2]),2); + l45=pow((x[4]-x[3]),2)+pow((y[4]-y[3]),2)+pow((z[4]-z[3]),2); + + cos_theta_triangle1=(l15+l12-l25)/(2*sqrt(l12*l15)); + cos_theta_triangle2=(l25+l23-l35)/(2*sqrt(l23*l25)); + cos_theta_triangle3=(l35+l34-l45)/(2*sqrt(l34*l35)); + cos_theta_triangle4=(l45+l14-l15)/(2*sqrt(l14*l45)); + + //First triangle : nodes 1, 2 and 5 + xyz_tria[0][0]=0; + xyz_tria[0][1]=0; + xyz_tria[1][0]=sqrt(l12); + xyz_tria[1][1]=0; + xyz_tria[2][0]=cos_theta_triangle1*sqrt(l15); + xyz_tria[2][1]=sqrt(1-pow(cos_theta_triangle1,2))*sqrt(l15); + + //Second triangle : nodes 2, 3 and 5 + xyz_tria[3][0]=0; + xyz_tria[3][1]=0; + xyz_tria[4][0]=sqrt(l23); + xyz_tria[4][1]=0; + xyz_tria[5][0]=cos_theta_triangle2*sqrt(l25); + xyz_tria[5][1]=sqrt(1-pow(cos_theta_triangle2,2))*sqrt(l25); + + //Third triangle : nodes 3, 4 and 5 + xyz_tria[6][0]=0; + xyz_tria[6][1]=0; + xyz_tria[7][0]=sqrt(l34); + xyz_tria[7][1]=0; + xyz_tria[8][0]=cos_theta_triangle3*sqrt(l35); + xyz_tria[8][1]=sqrt(1-pow(cos_theta_triangle3,2))*sqrt(l35); + + //Fourth triangle : nodes 4, 1 and 5 + xyz_tria[9][0]=0; + xyz_tria[9][1]=0; + xyz_tria[10][0]=sqrt(l14); + xyz_tria[10][1]=0; + xyz_tria[11][0]=cos_theta_triangle4*sqrt(l45); + xyz_tria[11][1]=sqrt(1-pow(cos_theta_triangle4,2))*sqrt(l45); + + + + //Start looping on the triangle's gaussian points: + for(ig=0;igterms[0]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*nx[i]; + pe_g->terms[1]+= J[i]*gauss_weight* pressure_tria[0]*l1l4_tria[i][0]*ny[i]; + pe_g->terms[2]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*nx[i]; + pe_g->terms[3]+= J[i]*gauss_weight* pressure_tria[1]*l1l4_tria[i][1]*ny[i]; + pe_g->terms[4]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*nx[i]; + pe_g->terms[5]+= J[i]*gauss_weight* pressure_tria[2]*l1l4_tria[i][2]*ny[i]; + pe_g->terms[6]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*nx[i]; + pe_g->terms[7]+= J[i]*gauss_weight* pressure_tria[3]*l1l4_tria[i][3]*ny[i]; + + + + } //for(i=0;i<4;i++) + } //for(ig=0;ig +#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" + +#include "./Matice.h" + + +/*-------------------------------------------------- + NewMatice + --------------------------------------------------*/ + +Matice* NewMatice(void) { + + return (Matice*)xmalloc(sizeof(Matice)); +} + + +/*-------------------------------------------------- + MaticeInit + --------------------------------------------------*/ + +int MaticeInit(Matice* this) +{ + + this->name[0] = '\0'; + this->mid = UNDEF; + this->B = UNDEF; + this->n = UNDEF; + + return 1; +} + + +/*-------------------------------------------------- + DeleteMatice + --------------------------------------------------*/ + +void DeleteMatice(void* *vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + Matice* *pthis = (Matice**)vpthis; + Matice* this = *pthis; + + #ifdef _DEBUG_ + _printf_(" DeleteMatice message:\n"); + _printf_(" deleting Matice with mid = %d\n", this->mid); + #endif + + xfree(vpthis); +} + + +/*-------------------------------------------------- + MaticeEcho + --------------------------------------------------*/ + +void MaticeEcho(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + Matice* this = (Matice*)vpthis; + + if (this) { + + printf("\nMatice record:\n"); + printf(" name: %s\n", this->name); + + if ( UNDEF!=this->mid) printf(" mid: %d\n", this->mid ); + else printf(" mid: UNDEF\n"); + + if( UNDEF!=(float)this->B) printf(" B: %lf\n",this->B ); + else printf(" B: UNDEF\n"); + + if( UNDEF!=(float)this->n)printf( " n: %g\n",this->n ); + else printf(" n: UNDEF\n"); + + } + else { + + printf(" null pointer input for Matice object; no echo\n"); + } +} + + +/*-------------------------------------------------- + MaticeSizeOfInDoubles + --------------------------------------------------*/ + +int MaticeSizeOfInDoubles(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + Matice* this = (Matice*)vpthis; + int struct_doubles=0; + + struct_doubles = + (int)ceil((float)sizeof(*this)/(float)sizeof(double)); + + #ifdef _DEBUG_ + _printf_("\n MaticeSizeOfInDoubles diagnostics:\n"); + _printf_(" returning %d doubles\n", struct_doubles); + #endif + + return struct_doubles; +} + + +/*-------------------------------------------------- + MaticeShrinkWrap + --------------------------------------------------*/ + +void MaticeShrinkWrap(void* to, void* vpthis) +{ + /* no "extended" data means a simple copy will do: */ + + memcpy( to, vpthis, sizeof(double)*MaticeSizeOfInDoubles(vpthis)); +} + + +/*-------------------------------------------------- + MaticeUnwrap + --------------------------------------------------*/ + +int MaticeUnwrap(void* vpthis) +{ + /* no extended data means no internal pointers to reestablish */ + + return 1; +} + + +/*-------------------------------------------------- + MaticeAddToDataSet + --------------------------------------------------*/ + +int MaticeAddToDataSet( DataSet* dataset, Matice* 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, MaticeEnum(), (void*)this, + &DeleteMatice, &MaticeEcho, &MaticeSizeOfInDoubles, + &MaticeShrinkWrap, &MaticeUnwrap, + &MaticeMarshall, + &MaticeGetID, &MaticeSetVirtualFunctions); + + return 1; +} + + +/*-------------------------------------------------- + MaticeGetID + --------------------------------------------------*/ + +int MaticeGetID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + + Matice* this = (Matice*)vpthis; + + return this->mid; +} + + +/*-------------------------------------------------- + MaticeGetName + --------------------------------------------------*/ + +char* MaticeGetName( void* vpthis) +{ + Matice* this = (Matice*)vpthis; + + #ifdef _DEBUG_ + _printf_("virtual function: MaticeGetName\n"); + #endif + + return this->name; +} + + +/*-------------------------------------------------- + MaticeSetVirtualFunctions + --------------------------------------------------*/ + +int MaticeSetVirtualFunctions( void* vp) +{ + + /* VirtualFunctions implemented as a "friend" class: */ + + /* void* for polymorphic function compatibility */ + VirtualFunctions* vf = (VirtualFunctions*)vp; + + if ( VirtualFunctionsInit(vf) ) { + + /* set Matice-specific virtual function pointers: */ + + /* general: */ + + vf->GetName = &MaticeGetName; + + return 1; + } + else + return 0; + +} + + +/*-------------------------------------------------- + MaticeMarshall + --------------------------------------------------*/ + +void MaticeMarshall( char* *pbuffer, void* vpthis, int* size) +{ + /* vpthis for polymorphic function compatibility */ + Matice* this = (Matice*)vpthis; + + char* buffer=*pbuffer; + + int size_string1; + + size_string1=strlen(this->name)+1; + + cellmar(&size_string1,&buffer,INTEGER_FLAG); + cellmar(&this->name,&buffer,STRING_FLAG); + cellmar(&this->mid,&buffer,INTEGER_FLAG); + cellmar(&this->B,&buffer,DOUBLE_FLAG); + cellmar(&this->n,&buffer,DOUBLE_FLAG); + + /* return values: */ + + *size=size_string1*sizeof(char)+ 2*sizeof(int)+2*sizeof(double); + *pbuffer=buffer; + + return; +} + + +/*-------------------------------------------------- + MaticeDemarshall + --------------------------------------------------*/ + +void* MaticeDemarshall( DataSet* ds, char* *pbuffer, int machine_flag) +{ + int size_string1; + + char* buffer=*pbuffer; + Matice* this=NewMatice(); + + celldem(&size_string1,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->name,&buffer,STRING_FLAG,machine_flag,size_string1); + celldem(&this->mid,&buffer,INTEGER_FLAG,machine_flag); + celldem(&this->B,&buffer,DOUBLE_FLAG,machine_flag); + celldem(&this->n,&buffer,DOUBLE_FLAG,machine_flag); + + if (ds) MaticeAddToDataSet( ds, this); + *pbuffer=buffer; + + return (void*)this; +} + +/*-------------------------------------------------- + MaticeGetViscosity2d + --------------------------------------------------*/ +int MaticeGetViscosity2d(double* pmu, void* vpthis, double* epsilon){ + +/*From a string tensor and a material object, return viscosity, using Glen's flow law. + 2*B + mu= ------------------------------------------------------------------- + 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + + where mu is the viscotiy, B the flow law parameter , (u,v) the velocity + vector, and n the flow law exponent. + + If epsilon is NULL, it means this is the first time Emg is being run, and we + return 10^14, initial viscosity. +*/ + + int noerr=1; + + /*output: */ + double mu; + + /*input from the material card: */ + Matice* this=NULL; /* vpthis for polymorphic function compatibility */ + double B,n; + + /*input strain rate: */ + double exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + double A,e; + + /*Retrieve material card info: */ + this=(Matice*)vpthis; + + + if (this->n==1){ + /*Viscous behaviour! mu=B: */ + mu=this->B; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){ + mu=pow(10,14); + } + else{ + B=this->B; + n=this->n; + + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: mu=2*B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity for 0 shear areas: */ + mu=4.5*pow(10,17); + } + else{ + e=(n-1)/2/n; + + mu=2*B/(2*pow(A,e)); + } + } + } + #ifdef _DEBUG_ + _printf_("Viscosity %lf\n",mu); + #endif + + /*Return: */ + *pmu=mu; + return noerr; +} + +/*-------------------------------------------------- + MaticeGetViscosity3d + --------------------------------------------------*/ +int MaticeGetViscosity3d(double* pmu, void* vpthis, double* epsilon){ + + /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: + * + * 2*B + * mu= ------------------------------------------------------------------- + * 2[ exx^2+eyy^2+exx*eyy+exy^2+exz^2+eyz^2 ]^[(n-1)/2n] + * + * where mu is the viscotiy, B the flow law parameter , (u,v) the velocity + * vector, and n the flow law exponent. + * + * If epsilon is NULL, it means this is the first time Emg is being run, and we + * return g, initial viscosity. + */ + + int noerr=1; + + /*output: */ + double mu; + + /*input from the material card: */ + Matice* this=NULL; /* vpthis for polymorphic function compatibility */ + double B,n; + + /*input strain rate: */ + double exx,eyy,exy,exz,eyz; + + /*Intermediary value A and exponent e: */ + double A,e; + + /*Retrieve material card info: */ + this=(Matice*)vpthis; + + if (this->n==1){ + /*Viscous behaviour! mu=B: */ + mu=this->B; + } + else{ + if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && + (epsilon[3]==0) && (epsilon[4]==0)){ + mu=pow(10,14); + } + else{ + + B=this->B; + n=this->n; + + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + exz=*(epsilon+3); + eyz=*(epsilon+4); + + /*Build viscosity: mu=2*B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+pow(exz,2)+pow(eyz,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity for 0 shear areas: */ + mu=4.5*pow(10,17); + } + else{ + e=(n-1)/2/n; + + mu=2*B/(2*pow(A,e)); + } + } + } + #ifdef _DEBUG_ + _printf_("Viscosity %lf\n",mu); + #endif + + /*Return: */ + *pmu=mu; + return noerr; +} + +/*-------------------------------------------------- + MaticeGetViscosity2 + --------------------------------------------------*/ +int MaticeGetViscosity2(double* pmu2, void* vpthis, double* epsilon){ + + /*Return viscosity accounting for steady state power law creep [Thomas and MacAyeal, 1982]: + * + * 2* (1-n)/2n + * mu2= ------------------------------------------------------------------- + * 2[ (du/dx)^2+(dv/dy)^2+1/4*(du/dy+dv/dx)^2+du/dx*dv/dy ]^[(3n-1)/2n] + * + * where mu2 is the second viscosity, (u,v) the velocity + * vector, and n the flow law exponent. + * + * If epsilon is NULL, it means this is the first time Gradjb is being run, and we + * return mu20, initial viscosity. + */ + + int noerr=1; + + /*output: */ + double mu2; + + /*input from the material card: */ + Matice* this=NULL; /* vpthis for polymorphic function compatibility */ + double n; + + /*input strain rate: */ + double exx,eyy,exy; + + /*Intermediary value A and exponent e: */ + double A,e; + + /*Retrieve material card info: */ + this=(Matice*)vpthis; + + if(epsilon){ + n=this->n; + + /*Retrive strain rate components: */ + exx=*(epsilon+0); + eyy=*(epsilon+1); + exy=*(epsilon+2); + + /*Build viscosity: mu2=B/(2*A^e) */ + A=pow(exx,2)+pow(eyy,2)+pow(exy,2)+exx*eyy; + if(A==0){ + /*Maximum viscosity2 for 0 shear areas: */ + mu2=4.5*pow(10,17); + } + else{ + e=(3*n-1)/2/n; + + mu2=2*(1-n)/2/n/(2*pow(A,e)); + } + } + else{ + mu2=4.5*pow(10,17); + } + + #ifdef _DEBUG_ + _printf_("Viscosity2 %lf\n",mu2); + #endif + + /*Return: */ + *pmu2=mu2; + return noerr; +} + +/*-------------------------------------------------- + MaticeSetFlowLawParameter + --------------------------------------------------*/ +int MaticeSetFlowLawParameter(void* vpthis, double B_value){ + + Matice* this=NULL; /* vpthis for polymorphic function compatibility */ + this=(Matice*)vpthis; + + this->B=B_value; + return 1; +} Index: /issm/trunk/src/c/objects/cielo/Matice.h =================================================================== --- /issm/trunk/src/c/objects/cielo/Matice.h (revision 216) +++ /issm/trunk/src/c/objects/cielo/Matice.h (revision 216) @@ -0,0 +1,60 @@ +/* + * Matice.h + */ + + +#ifndef _MATICE_H +#define _MATICE_H + + +/* "friend" functionality: */ +#include "../DataSet/DataSet.h" + +/* input file-specific defines: */ +#include "./InputFormats.h" + +typedef struct { + char name[LARGEFIELDWIDTH+1]; + int mid; + double B; + double n; +} Matice; + + +/* creation, initialisation: */ + + Matice* NewMatice(void); + int MaticeInit(Matice* this); + + +/* "virtual functions", used or required by DataSet: */ + + void DeleteMatice( void* *vpthis); /* void*'s only for compatibility */ + void MaticeEcho( void* vpthis); + int MaticeSizeOfInDoubles( void* vpthis); + void MaticeShrinkWrap( void* to, void* vpthis); + int MaticeUnwrap( void* vpthis); + void MaticeMarshall( char* *pbuffer, void* vpthis, int* size); + void* MaticeDemarshall( DataSet* ds, char* *pbuffer, int machine_flag); + int MaticeGetID( void* this); + int MaticeSetVirtualFunctions( void* vf); + + +/* other member functions: */ + + /* alternative interface to DataSetAddObject */ + int MaticeAddToDataSet( DataSet* dataset, Matice* this); + int MaticeGetViscosity2d(double* pmu, void* vpthis, double* epsilon); + int MaticeGetViscosity3d(double* pmu, void* vpthis, double* epsilon); + int MaticeSetFlowLawParameter(void* vpthis,double B_value); + int MaticeGetViscosity2(double* pmu2, void* vpthis, double* epsilon); + +/* other "virtual" functions: */ + + /* general: */ + + char* MaticeGetName( void* vpthis); + + +#endif /* _MATICE_H */ + Index: /issm/trunk/src/c/objects/cielo/PengridLoad.c =================================================================== --- /issm/trunk/src/c/objects/cielo/PengridLoad.c (revision 216) +++ /issm/trunk/src/c/objects/cielo/PengridLoad.c (revision 216) @@ -0,0 +1,1157 @@ + +/* + PengridLoad.c: +*/ + +#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 "./PengridLoad.h" + +/*Objects: */ +#include "../objects.h" + +/*-------------------------------------------------- + NewPengridLoad + --------------------------------------------------*/ + +PengridLoad* NewPengridLoad( void* vppengrid) +{ + + Pengrid* pengrid = (Pengrid*)vppengrid; + PengridLoad* this = NULL; + int i,j; + + this=(PengridLoad*)xcalloc(1,sizeof(PengridLoad)); /* guaranteed NULL, 0 initialisation */ + + /* include the "base class" in the "derived" object: */ + memcpy( &this->pengrid, pengrid, sizeof(Pengrid)); + + /* since zero could be a valid index, set all DataSet index + "pointers" to UNDEF: */ + this->internalgrid_index=UNDEF; + this->internalgrid_indexb=UNDEF; + this->pg=NULL; + + this->matpar = NULL; + this->matpar_index=0; + this->matpar_enum=0; + + return this; +} + + +/*-------------------------------------------------- + DeletePengridLoad + --------------------------------------------------*/ + +void DeletePengridLoad(void* *vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + PengridLoad* *pthis = (PengridLoad**)vpthis; + PengridLoad* this = *pthis; + + xfree(vpthis); +} + + +/*-------------------------------------------------- + PengridLoadEcho + --------------------------------------------------*/ + +void PengridLoadEcho(void* vpthis){ + int i; + + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = (PengridLoad*)vpthis; + + if (this) { + + _printf_("\nPengridLoad echo:\n"); + printf("-------------------\n"); + printf(" base Pengrid object:\n"); + PengridEcho( &this->pengrid); + printf("\n"); + _printf_(" penalty grid points to following grids (pointer-based ops):\n"); + _printf_(" local grid number corresponding grid id\n"); + _printf_(" ----------------- ---------------------\n"); + if ( this->pg) + _printf_(" %8d %8d\n", i, InternalGridGetID( this->pg)); + if ( this->matpar) printf(" matpar entry of enum type %d\n", this->matpar_enum); + else printf(" matpar reference is NULL\n"); + + } + else{ + _printf_(" null pointer input for PengridLoad object; no echo\n"); + } +} + +/*-------------------------------------------------- + PengridLoadSizeOfInDoubles + --------------------------------------------------*/ + +int PengridLoadSizeOfInDoubles(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + PengridLoad* this = (PengridLoad*)vpthis; + int struct_doubles=0; + + struct_doubles = + (int)ceil((float)sizeof(*this)/(float)sizeof(double)); + + #ifdef _DEBUG_ + _printf_("\n PengridLoadSizeOfInDoubles diagnostics:\n"); + _printf_(" returning %d doubles\n", struct_doubles); + #endif + + return struct_doubles; +} + + +/*-------------------------------------------------- + PengridLoadShrinkWrap + --------------------------------------------------*/ + +void PengridLoadShrinkWrap(void* to, void* vpthis) +{ + /* no "extended" data means a simple copy will do: */ + + memcpy( to, vpthis, sizeof(double)*PengridLoadSizeOfInDoubles(vpthis)); +} + + +/*-------------------------------------------------- + PengridLoadUnwrap + --------------------------------------------------*/ + +int PengridLoadUnwrap(void* vpthis) +{ + PengridLoad* this = (PengridLoad*)vpthis; + + this->matpar = NULL; + + return 1; +} + +/*-------------------------------------------------- + PengridLoadMarshall + --------------------------------------------------*/ + +void PengridLoadMarshall( char* *pbuffer, void* vpthis, int* size) +{ + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = (PengridLoad*)vpthis; + + char* buffer=*pbuffer; + + int i, j, size_pengrid=0; + int ig; + + int num_ints=0, num_doubles=0; + + + /* first marshall the embedded Pengrid object: ... */ + + PengridMarshall( &buffer, &this->pengrid, &size_pengrid); + + + /* ... then marshall the PengridLoad member data: */ + + cellmar( &this->internalgrid_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->internalgrid_indexb, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->pg, &buffer, INTEGER_FLAG); num_ints++; + + cellmar( &this->matpar_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matpar_enum, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matpar, &buffer, INTEGER_FLAG); num_ints++; + + *size = size_pengrid + num_ints*sizeof(int) + num_doubles*sizeof(double); + *pbuffer = buffer; + + return; +} + + +/*-------------------------------------------------- + PengridLoadDemarshall + --------------------------------------------------*/ + +void* PengridLoadDemarshall( DataSet* ds, char* *pbuffer, int machine_flag) +{ + int i, j, ig; + Pengrid* ppengrid=NULL; /* for embedded object */ + + char* buffer=*pbuffer; + PengridLoad* this=NULL; + + + /* first demarshall the embedded Pengrid object: ... */ + + ppengrid = PengridDemarshall( NULL, &buffer, machine_flag); + /* create the derived object: */ + this = NewPengridLoad( ppengrid); + /* and free pengrid now that it's embedded in the derived object: */ + DeletePengrid( (void**)&ppengrid); + + + /* ... then demarshall the PengridLoad member data: */ + + celldem( &this->internalgrid_index, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->internalgrid_indexb, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->pg, &buffer, INTEGER_FLAG, machine_flag); + + celldem( &this->matpar_index, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->matpar_enum, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->matpar, &buffer, INTEGER_FLAG, machine_flag); + + /*Set pointers to NULL :*/ + this->pg=NULL; + this->matpar=NULL; + + if (ds) PengridLoadAddToDataSet( ds, this); + *pbuffer=buffer; + + return (void*)this; +} + +/*-------------------------------------------------- + PengridLoadGetID + --------------------------------------------------*/ + +int PengridLoadGetID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = (PengridLoad*)vpthis; + + return PengridGetID( &this->pengrid); +} + +/*-------------------------------------------------- + PengridLoadSetVirtualFunctions + --------------------------------------------------*/ + +int PengridLoadSetVirtualFunctions( void* vp) +{ + + /* VirtualFunctions implemented as a "friend" class: */ + + /* void* for polymorphic function compatibility */ + VirtualFunctions* vf = (VirtualFunctions*)vp; + + if ( VirtualFunctionsInit(vf) ) { + + /* set PengridLoad-specific virtual function pointers: */ + + /* general: */ + + /* configuration: */ + vf->Configure = &PengridLoadConfigure; + + /* penalty stiffnesses and loads: */ + vf->CreateKMatrix = NULL; + vf->PenaltyCreateKMatrix = &PengridLoadPenaltyCreateKMatrix; + vf->CreatePVector = NULL; + vf->PenaltyCreatePVector = &PengridLoadPenaltyCreatePVector; + vf->PenaltyConstrain = &PengridLoadPenaltyConstrain; + + return 1; + } + else + return 0; + +} + +/*-------------------------------------------------- + PengridLoadAddToDataSet + --------------------------------------------------*/ + +int PengridLoadAddToDataSet( DataSet* dataset, PengridLoad* 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, PengridLoadEnum(), (void*)this, + &DeletePengridLoad, &PengridLoadEcho, &PengridLoadSizeOfInDoubles, + &PengridLoadShrinkWrap, &PengridLoadUnwrap, + &PengridLoadMarshall, + &PengridLoadGetID, &PengridLoadSetVirtualFunctions); + + return 1; +} + + +/*-------------------------------------------------- + PengridLoadConfigure + --------------------------------------------------*/ + +int PengridLoadConfigure( void* vpthis, int num_datasets, ...) { + + /* set up data connections to all necessary, related objects: */ + + PengridLoad* this = (PengridLoad*)vpthis; + + int i,j, found, foundgrid, noerr=1, first=1; + + /* expected DataSet input: */ + + DataSet* pgpdt=NULL; + DataSet* pgpdtb=NULL; + DataSet* pmpt=NULL; + + /* variable argument list handling: */ + + va_list args; + DataSet* *arglist=NULL; + + #ifdef _DEBUG_ + _printf_("virtual function: PengridLoadConfigure\n"); + #endif + if (this) { + + if (!num_datasets) { + + _printf_(" %s error:\n", __func__); + _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*/ + this->pg=NULL; + + /* grid connections based on DataSet of type "gpdt": */ + + #ifndef _PARALLEL_ + found=0; + for ( i=0; ipengrid); + + if ( (current_grid_id=pgrid_id_list[0]) ) { /* 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_index ) { + + /* does index still point to a valid grid? (i.e., the one with + the same id): */ + + if ( InternalGridEnum() == DataSetGetEnum( pgpdt, this->internalgrid_index) + && + current_grid_id == DataSetGetObjectID( pgpdt, this->internalgrid_index) ) { + + /* great, still valid! since we have the index position in + gpdt, just use it to get, save, the grid pointer: */ + + this->pg = DataSetGetObjectPtr( pgpdt, this->internalgrid_index); + } + } + + /* if we haven't yet "recovered" the InternalGrid object pointer, + go seek it out: */ + + if ( !this->pg ) { + + noerr *= DataSetSeek( (void**)(&this->pg), /* cast because some compilers recognise as */ + /* ptr to ptr to struct otherwise */ + &this->internalgrid_index , + pgpdt, current_grid_id, 1, InternalGridEnum()); + + if (!noerr) { + if (first) { + _printf_(" %s error:\n", __func__); + first = 0; + } + _printf_(" could not find grid with id of %d for load %d of type \"%s\"\n", + current_grid_id, PengridLoadGetID(this), PengridLoadGetName(this)); + } + } + + } + } + else { + + if (first) { + _printf_(" %s error:\n", __func__); + first=0; + } + _printf_(" could not find DataSet of type \"gpdt\", grid object\n"); + _printf_(" pointers not established for pengrid load id = %d\n", + PengridLoadGetID(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; ipengrid); + + if ( (current_grid_id=pgrid_id_list[0]) ) { /* 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_index ) { + + /* does index still point to a valid grid? (i.e., the one with + the same id): */ + + if ( InternalGridEnum() == DataSetGetEnum( pgpdt, this->internalgrid_index) + && + current_grid_id == DataSetGetObjectID( pgpdt, this->internalgrid_index) ) { + + /* great, still valid! since we have the index position in + gpdt, just use it to get, save, the grid pointer: */ + + this->pg = DataSetGetObjectPtr( pgpdt, this->internalgrid_index); + } + } + if ( UNDEF != (int)this->internalgrid_indexb ) { + + /* does index still point to a valid grid? (i.e., the one with + the same id): */ + + if ( InternalGridEnum() == DataSetGetEnum( pgpdtb, this->internalgrid_indexb) + && + current_grid_id == DataSetGetObjectID( pgpdtb, this->internalgrid_indexb) ) { + + /* great, still valid! since we have the index position in + gpdtb, just use it to get, save, the grid pointer: */ + + this->pg = DataSetGetObjectPtr( pgpdtb, this->internalgrid_indexb); + } + } + + /* if we haven't yet "recovered" the InternalGrid object pointer, + go seek it out: */ + + if ( !this->pg ) { + + foundgrid = DataSetSeek( (void**)(&this->pg), /* cast because some compilers recognise as */ + /* ptr to ptr to struct otherwise */ + &this->internalgrid_index , + 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), /* cast because some compilers recognise as */ + /* ptr to ptr to struct otherwise */ + &this->internalgrid_indexb , + 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, PengridLoadGetID(this), PengridLoadGetName(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 pengrid element id = %d\n", + PengridLoadGetID(this)); + noerr = 0; + } + + #endif + + + + + + /* matice property connections based on DataSet of type "mpt": */ + found=0; + for ( i=0; ipengrid); + + /* as for grids, check to see if a nonzero matpar index still + points to the matpar entry of interest in mpt, otherwise do + a seek for it: */ + + if ( this->matpar_enum != UNDEF ){ + + + if ( this->matpar_enum == DataSetGetEnum( pmpt, this->matpar_index) + && + current_matpar_id == DataSetGetObjectID( pmpt, this->matpar_index) ) { + + /* index is still valid; just get matpar entry pointer directly: */ + + this->matpar = DataSetGetObjectPtr( pmpt, this->matpar_index); + } + } + + /* if we haven't yet "recovered" the matpar object pointer, go seek it out: */ + + if ( NULL == this->matpar) { + + noerr *= DataSetSeek( &this->matpar, &this->matpar_index, + pmpt, current_matpar_id, 1, MatparEnum() ); /* 0 to generalise later to Pcomp, etc. ... */ + + if (!noerr) { + if (first) { + _printf_(" %s error:\n", __func__); + first = 0; + } + _printf_(" could not find matpar entry with id of %d for element %d of type \"%s\"\n", + current_matpar_id, PengridGetID(&this->pengrid), PengridGetName(&this->pengrid)); + } + else { + + /* though the answer for now is just MaticeEnum(), for generality go ahead and + ask the question, "what type of matpar has been referenced?" */ + + this->matpar_enum = DataSetGetEnum( *(arglist+i), this->matpar_index); + } + } + } + else { + + if (first) { + _printf_(" %s error:\n", __func__); + first=0; + } + _printf_(" could not find DataSet of type \"mpt\"; matpar relationships\n"); + _printf_(" not established for penta element id = %d\n", + PentaElementGetID(this)); + noerr = 0; + } + + #ifdef _DEBUG_ + _printf_(" load data after all lookups:\n"); + PengridLoadEcho( this); + #endif + + xfree((void**)&arglist); + + } + else { + + /* quiet return */ + ; + } + cleanup_and_return: + + return noerr; + +} +/*-------------------------------------------------- + PengridLoadPenaltyCreateKMatrix + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyCreateKMatrix" +int PengridLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){ + + int noerr=1; + if ((analysis_type==ThermalSteadyAnalysisEnum()) || (analysis_type==ThermalTransientAnalysisEnum())){ + + PengridLoadPenaltyCreateKMatrixThermal( pKe_gg, vpthis, inputs, K_flag, kmax, analysis_type); + } + else if (analysis_type==DiagnosticStokesAnalysisEnum()){ + _printf_("%s%s\n",__FUNCT__," error message: analysis_type not supported yet!\n"); + return 0; + } + else if (analysis_type==MeltingAnalysisEnum()){ + + PengridLoadPenaltyCreateKMatrixMelting( pKe_gg, vpthis, inputs, K_flag, kmax, analysis_type); + } + else{ + _printf_("%s%s\n",__FUNCT__," error message: analysis_type not supported yet!\n"); + return 0; + } + return noerr; + +} +/*-------------------------------------------------- + PengridLoadPenaltyCreateKMatrixThermal + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyCreateKMatrixThermal" +int PengridLoadPenaltyCreateKMatrixThermal( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){ + + int i,counter; + int noerr=1; + + PengridLoad* this = NULL; + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int* structural_dof_list = NULL; + int dof; + + double* dt_param; + double dt; + + + /*Some pointer intialization: */ + this = (PengridLoad*)vpthis; + + /*Is penalty active? :*/ + if (this->pengrid.active){ + + /*Recover dt input: */ + dt_param=ParameterInputsRecover(inputs,"dt"); + if (!dt_param){ + _printf_("%s%s\n",__FUNCT__," error message: dt needed to apply penalty!"); + *pKe_gg=NULL; + return 0; + } + else{ + dt=*dt_param; + } + + /*initialize stiffness matrix: */ + Ke_gg=NewElemMatrix(1); + + /*Build dof list for stiffness matrix: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg); + dof=structural_dof_list[this->pengrid.dof-1]; + Ke_gg->row_indices[0]=dof; + + if (this->pengrid.thermal_steadystate){ + *(Ke_gg->terms+Ke_gg->nrows*0+0)=this->pengrid.penalty_melting; + } + else{ + *(Ke_gg->terms+Ke_gg->nrows*0+0)=dt*this->pengrid.penalty_melting; + } + + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + } + else{ + Ke_gg=NULL; + } + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + return noerr; +} +/*-------------------------------------------------- + PengridLoadPenaltyCreateKMatrixMelting + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyCreateKMatrixMelting" +int PengridLoadPenaltyCreateKMatrixMelting( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type){ + + int i,counter; + int noerr=1; + + PengridLoad* this = NULL; + Matpar* matpar = NULL; + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int* structural_dof_list = NULL; + int dof; + + double* temperature_param; + double temperature; + + double* pressure_param; + double pressure; + double meltingpoint,beta,t_pmp; + + /*Some pointer intialization: */ + this = (PengridLoad*)vpthis; + matpar=this->matpar; + + + //Recover parameters: */ + meltingpoint=matpar->meltingpoint; + beta=matpar->beta; + + //First recover pressure and temperature vectors from inputs + temperature_param=ParameterInputsRecover(inputs,"temperature"); + pressure_param=ParameterInputsRecover(inputs,"pressure"); + + if((!pressure_param) || (!temperature_param)){ + _printf_("%s%s\n",__FUNCT__," error message: missing input temperature or pressure!"); + return 0; + } + + /*Recover pressure and temperatureat the grid:*/ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg); + dof=structural_dof_list[0]; + temperature=temperature_param[dof]; + pressure=pressure_param[dof]; + + /*Compute pressure melting point: */ + t_pmp=meltingpoint-beta*pressure; + + + /*Add penalty load: */ + if (temperature>t_pmp){ + Ke_gg=NULL; + } + else{ + //If Trow_indices[0]=dof; + *(Ke_gg->terms+Ke_gg->nrows*0+0)=kmax*pow(10,this->pengrid.penalty_offset); + } + + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return noerr; +} + +/*-------------------------------------------------- + PengridLoadPenaltyCreatePVector + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyCreatePVector" + +int PengridLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ + + int noerr=1; + if ((analysis_type==ThermalSteadyAnalysisEnum()) || (analysis_type==ThermalTransientAnalysisEnum())){ + + PengridLoadPenaltyCreatePVectorThermal(ppe_g, vpthis, inputs, kmax, analysis_type); + } + else if (analysis_type==DiagnosticStokesAnalysisEnum()){ + _printf_("%s%s\n",__FUNCT__," error message: analysis_type not supported yet!\n"); + return 0; + } + else if (analysis_type==MeltingAnalysisEnum()){ + + PengridLoadPenaltyCreatePVectorMelting(ppe_g, vpthis, inputs, kmax, analysis_type); + } + else{ + _printf_("%s%s\n",__FUNCT__," error message: analysis_type not supported yet!\n"); + return 0; + } + return noerr; +} + + +/*-------------------------------------------------- + PengridLoadPenaltyCreatePVectorThermal + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyCreatePVectorThermal" + +int PengridLoadPenaltyCreatePVectorThermal( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ + + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = NULL; + Matpar* matpar = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + /*input parameters: */ + double* dt_param=NULL; + double dt; + double* pressure_param=NULL; + double pressure; + + /*material parameters: */ + double beta; + double meltingpoint; + double t_pmp; + + int* structural_dof_list = NULL; + int dof; + + /*Some pointer intialization: */ + this = (PengridLoad*)vpthis; + matpar=this->matpar; + + /*Is penalty load active? :*/ + if (!this->pengrid.active){ + /*do nothing: */ + pe_g=NULL; + } + else{ + + /* Allocate new load vector: */ + pe_g=NewElemVector(1); + + /* recover input parameters: */ + dt_param=ParameterInputsRecover(inputs,"dt"); + pressure_param=ParameterInputsRecover(inputs,"pressure"); + + if( (!pressure_param)){ + _printf_("%s%s\n",__FUNCT__," error message: missing pressure in input parameters!"); + noerr=0;goto cleanup_and_return; + } + if(this->pengrid.thermal_steadystate==0){ + if( (!dt_param)){ + _printf_("%s%s\n",__FUNCT__," error message: missing dt in input parameters!"); + noerr=0;goto cleanup_and_return; + } + } + + /*Recover pressure at the grid: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr( this->pg); + dof=structural_dof_list[0]; //thermal dof #1 + *(pe_g->row_indices+0)=dof; + pressure=pressure_param[dof]; + + /*Recover material parameters: */ + beta=matpar->beta; + meltingpoint=matpar->meltingpoint; + + //compute pressure melting point + t_pmp=meltingpoint-beta*pressure; + + //add penalty load: + if(this->pengrid.thermal_steadystate){ + *(pe_g->terms+0)= this->pengrid.penalty_melting*t_pmp; + } + else{ + dt=*dt_param; + *(pe_g->terms+0)= dt*this->pengrid.penalty_melting*t_pmp; + } + + + } + + cleanup_and_return: + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } + +} + +/*-------------------------------------------------- + PengridLoadPenaltyCreatePVectorMelting + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyCreatePVectorMelting" + +int PengridLoadPenaltyCreatePVectorMelting( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type){ + + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = NULL; + Matpar* matpar = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + /*input parameters: */ + double* pressure_param=NULL; + double pressure; + double* temperature_param=NULL; + double temperature; + + /*material parameters: */ + double beta; + double meltingpoint; + double heatcapacity; + double latentheat; + double t_pmp; + + int* structural_dof_list = NULL; + int dof; + + /*Some pointer intialization: */ + this = (PengridLoad*)vpthis; + matpar=this->matpar; + + /* Allocate new load vector: */ + pe_g=NewElemVector(1); + + /* recover input parameters: */ + pressure_param=ParameterInputsRecover(inputs,"pressure"); + temperature_param=ParameterInputsRecover(inputs,"temperature"); + + if( (!pressure_param) || (!temperature_param)){ + _printf_("%s%s\n",__FUNCT__," error message: missing pressure or temperature in input parameters!"); + noerr=0;goto cleanup_and_return; + } + + /*Recover pressure and temperature at the grid: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr( this->pg); + dof=structural_dof_list[0]; //thermal dof #1 + *(pe_g->row_indices+0)=dof; + pressure=pressure_param[dof]; + temperature=temperature_param[dof]; + + /*Recover material parameters: */ + beta=matpar->beta; + meltingpoint=matpar->meltingpoint; + heatcapacity=matpar->heatcapacity; + latentheat=matpar->latentheat; + + //compute pressure melting point + t_pmp=meltingpoint-beta*pressure; + + //add penalty load: + if (temperature>t_pmp){ + *(pe_g->terms+0)= this->pengrid.penalty_melting*(temperature-t_pmp); + } + else{ + *(pe_g->terms+0)= 0; + } + + cleanup_and_return: + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } + +} +/*-------------------------------------------------- + PengridLoadPenaltyConstrain + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadPenaltyConstrain" + +int PengridLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type){ + + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = NULL; + Matpar* matpar = NULL; + + int noerr=1; + + /*parameters: */ + double* pressure_param=NULL; + double pressure; + double* temperature_param=NULL; + double temperature; + + /*Material parameters: */ + double beta; + double meltingpoint; + double t_pmp; + + int* structural_dof_list = NULL; + int dof; + int unstable; + int active; + + /*Some pointer intialization: */ + this = (PengridLoad*)vpthis; + matpar=this->matpar; + + /*First recover parameter inputs: */ + pressure_param=ParameterInputsRecover(inputs,"pressure"); + temperature_param=ParameterInputsRecover(inputs,"temperature"); + + if(!pressure_param || !temperature_param ){ + _printf_("%s%s\n",__FUNCT__," error message: missing input parameters!"); + noerr=0 ;goto cleanup_and_return; + } + + /*Recover pressure and temperature at grid: */ + structural_dof_list= InternalGridGetStructuralDoflistPtr(this->pg); + dof=structural_dof_list[0]; //that's where temperatures and pressure were solved for. + temperature=temperature_param[dof]; + pressure=pressure_param[dof]; + + //Compute pressure melting point + meltingpoint=matpar->meltingpoint; + beta=matpar->beta; + t_pmp=meltingpoint-beta*pressure; + + //Figure out if temperature is over meltingpoint, in which case, this penalty needs to be activated. + if (temperature>t_pmp){ + active=1; + } + else{ + active=0; + } + + //Optional! + /*if(this->pengrid.active){ //once activated, a constraint will remain on + active=1; + }*/ + + //Figure out stability of this penalty + if(this->pengrid.active==active){ + unstable=0; + } + else{ + unstable=1; + } + + //Set penalty flag + this->pengrid.active=active; + + cleanup_and_return: + /*assign output pointer: */ + *punstable=unstable; + + return noerr; +} + + +/*-------------------------------------------------- + PengridLoadGetName + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PengridLoadGetName" + +char* PengridLoadGetName( void* vpthis){ + + /* vpthis for polymorphic function compatibility */ + PengridLoad* this = (PengridLoad*)vpthis; + + return PengridGetName(&this->pengrid); + +} Index: /issm/trunk/src/c/objects/cielo/PengridLoad.h =================================================================== --- /issm/trunk/src/c/objects/cielo/PengridLoad.h (revision 216) +++ /issm/trunk/src/c/objects/cielo/PengridLoad.h (revision 216) @@ -0,0 +1,74 @@ + +/* + PengridLoad.h +*/ + + +#ifndef _PENGRIDLOAD_H +#define _PENGRIDLOAD_H + + +#include "./PengridLoad.h" + +#include "../InternalGrid.h" + +/* element matrices: */ +#include "../objects.h" + +#include "../Pengrid.h" + +typedef struct { + + Pengrid pengrid; + + /* "friend" class object pointers and indices: */ + int internalgrid_index; + int internalgrid_indexb; + InternalGrid* pg; + + int matpar_index; + int matpar_enum; + void* matpar; + +} PengridLoad; + + +/* creation, initialisation: */ + + PengridLoad* NewPengridLoad( void* vppengrid); /* cast to Pengrid */ + + /* "virtual" functions, used or required by DataSet: */ + + void DeletePengridLoad( void* *vpthis); /* void*'s only for compatibility */ + void PengridLoadEcho( void* vpthis ); + int PengridLoadSizeOfInDoubles( void* vpthis); + void PengridLoadShrinkWrap( void* to, void* vpthis); + int PengridLoadUnwrap( void* vpthis); + void PengridLoadMarshall( char* *pbuffer, void* vpthis, int* size); + void* PengridLoadDemarshall( DataSet* ds, char* *pbuffer, int machine_flag); + int PengridLoadGetID( void* this); + int PengridLoadSetVirtualFunctions( void* vf); + + +/* other "virtual" functions: */ + + /* general: */ + + char* PengridLoadGetName( void* vpthis); + int PengridLoadConfigure( void* vpthis, int num_datasets, ...); + + int PengridLoadAddToDataSet( DataSet* dataset, PengridLoad* this); + + /* load functionality: */ + int PengridLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax,int analysis_type); + int PengridLoadPenaltyCreateKMatrixThermal( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type); + int PengridLoadPenaltyCreateKMatrixMelting( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type); + + int PengridLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type); + int PengridLoadPenaltyCreatePVectorThermal( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type); + int PengridLoadPenaltyCreatePVectorMelting( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type); + + int PengridLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type); + +#endif + Index: /issm/trunk/src/c/objects/cielo/PenpairLoad.h =================================================================== --- /issm/trunk/src/c/objects/cielo/PenpairLoad.h (revision 216) +++ /issm/trunk/src/c/objects/cielo/PenpairLoad.h (revision 216) @@ -0,0 +1,87 @@ + +/* + PenpairLoad.h +*/ + +/* + PenpairLoad is derived from the basic Penpair, which holds one penalty + pair, for 2 grids, for a certain dof. PenpairLoad adds load-level functionality + to the object. + + See PenpairLoad.c for more explanations. +*/ + + +#ifndef _PENPAIRLOAD_H +#define _PENPAIRLOAD_H + + +#include "./PenpairLoad.h" + +#include "../InternalGrid.h" + +/* element matrices: */ +#include "../objects.h" + + + + +typedef struct { + + Penpair penpair; + + /* "friend" class object pointers and indices: */ + int internalgrid_indices[2]; /* DataSet indices into partition grids dataset... */ + int internalgrid_indicesb[2]; /* DataSet indices into boundary partition grids dataset... */ + InternalGrid* pg[2]; /* ... and object pointers */ + + int active; + int counter; + int prestable; + + int internalelement_indices[2]; /* points to elements holding this penalty pair.*/ + void* elements[2]; + + +} PenpairLoad; + + +/* creation, initialisation: */ + + PenpairLoad* NewPenpairLoad( void* vppenpair); /* cast to Penpair */ + + /* "virtual" functions, used or required by DataSet: */ + + void DeletePenpairLoad( void* *vpthis); /* void*'s only for compatibility */ + void PenpairLoadEcho( void* vpthis ); + int PenpairLoadSizeOfInDoubles( void* vpthis); + void PenpairLoadShrinkWrap( void* to, void* vpthis); + int PenpairLoadUnwrap( void* vpthis); + void PenpairLoadMarshall( char* *pbuffer, void* vpthis, int* size); + void* PenpairLoadDemarshall( DataSet* ds, char* *pbuffer, int machine_flag); + int PenpairLoadGetID( void* this); + int PenpairLoadSetVirtualFunctions( void* vf); + + +/* other "virtual" functions: */ + + /* general: */ + + char* PenpairLoadGetName( void* vpthis); + int PenpairLoadConfigure( void* vpthis, int num_datasets, ...); + + int PenpairLoadAddToDataSet( DataSet* dataset, PenpairLoad* this); + + /* load functionality: */ + int PenpairLoadPenaltyCreateKMatrix( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type); + int PenpairLoadPenaltyCreateKMatrixOneDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inptus, int K_flag, double kmax, int analysis_type); + int PenpairLoadPenaltyCreateKMatrixTwoDof( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, double kmax, int analysis_type); + int PenpairLoadPenaltyConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type); + int PenpairLoadPenaltyCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type); + int PenpairLoadPenaltyCreatePVectorTwoDof( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, double kmax, int analysis_type); + int PenpairLoadPenetration(double* ppenetration, void* vpthis,ParameterInputs* inputs, int analysis_type); + int PenpairLoadPenaltyPreConstrain(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type); + int PenpairLoadPotentialUnstableConstraint(int* punstable, void* vpthis,ParameterInputs* inputs, int analysis_type); + +#endif + Index: /issm/trunk/src/c/objects/cielo/PentaElement.c =================================================================== --- /issm/trunk/src/c/objects/cielo/PentaElement.c (revision 216) +++ /issm/trunk/src/c/objects/cielo/PentaElement.c (revision 216) @@ -0,0 +1,5309 @@ +/* + PentaElement.c +*/ + + +#include +#include +#include + +/* variadic functions: */ +#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" + +#include "./PentaElement.h" + +#include "./GaussPoints.h" +#include "./GetElementGridData.h" +#include "./MatrixOperations.h" +#include "../Matice.h" + +/*Objects: */ +#include "../objects.h" + +#include "./PentaElement.h" + + +/*For debugging purposes: */ +#define ELID 141 //element for which to print debug statements +#define RANK 2 //rank of cpu for which to print debug statements. +//#define _DEBUGELEMENTS_ +//#define _DEBUGGAUSS_ + + +#define NDOF1 1 +#define NDOF2 2 +#define NDOF3 3 +#define numgrids 6 + +/*-------------------------------------------------- + NewPentaElement + --------------------------------------------------*/ +PentaElement* NewPentaElement( void* vppenta) +{ + + Penta* penta = (Penta*)vppenta; + PentaElement* this = NULL; + int i,j; + + this=(PentaElement*)xcalloc(1,sizeof(PentaElement)); /* guaranteed NULL, 0 initialisation */ + + /* include the "base class" in the "derived" object: */ + memcpy( &this->penta, penta, sizeof(Penta)); + + /* since zero could be a valid index, set all DataSet index + "pointers" to UNDEF: */ + + for ( i=0; iinternalgrid_indices[i] = UNDEF; + this->internalgrid_indicesb[i] = UNDEF; + this->pg[i]=NULL; + } + + this->matice = NULL; + this->matice_index=0; + this->matice_enum=0; + + this->matpar = NULL; + this->matpar_index=0; + this->matpar_enum=0; + + + return this; +} + + +/*-------------------------------------------------- + DeletePentaElement + --------------------------------------------------*/ + +void DeletePentaElement(void* *vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + PentaElement* *pthis = (PentaElement**)vpthis; + PentaElement* this = *pthis; + + #ifdef _DEBUG_ + _printf_(" DeletePentaElement message:\n"); + _printf_(" deleting PentaElement with eid = %d\n", this->penta.eid); + #endif + + xfree(vpthis); +} + +/*-------------------------------------------------- + PentaElementEcho + --------------------------------------------------*/ + +void PentaElementEcho(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + PentaElement* this = (PentaElement*)vpthis; + + int i,j; + ; + + printf("\nPentaElement echo:\n"); + printf("-------------------\n"); + printf(" base Penta object:\n"); + PentaEcho( &this->penta); + printf("\n"); + + printf(" element points to following grids (pointer-based ops):\n"); + + printf(" local grid number corresponding grid id\n"); + printf(" ----------------- ---------------------\n"); + + for ( i=0; ipg[i]) + printf(" %8d %8d\n", i, InternalGridGetID( this->pg[i])); + } + + printf(" element references the following:\n"); + + if ( this->matice) printf(" matice entry of enum type %d\n", this->matice_enum); + else printf(" matice reference is NULL\n"); + + if ( this->matpar) printf(" matpar entry of enum type %d\n", this->matpar_enum); + else printf(" matpar reference is NULL\n"); + +} + + +/*-------------------------------------------------- + PentaElementSizeOfInDoubles + --------------------------------------------------*/ + +int PentaElementSizeOfInDoubles(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + PentaElement* this = (PentaElement*)vpthis; + int struct_doubles=0; + + struct_doubles = + (int)ceil((float)sizeof(*this)/(float)sizeof(double)); + + #ifdef _DEBUG_ + _printf_("\n PentaElementSizeOfInDoubles diagnostics:\n"); + _printf_(" returning %d doubles\n", struct_doubles); + #endif + + return struct_doubles; +} + +/*-------------------------------------------------- + PentaElementShrinkWrap + --------------------------------------------------*/ + +void PentaElementShrinkWrap(void* to, void* vpthis) +{ + /* no "extended" data means a simple copy will do: */ + + memcpy( to, vpthis, PentaElementSizeOfInDoubles(vpthis)*sizeof(double)); +} + + +/*-------------------------------------------------- + PentaElementUnwrap + --------------------------------------------------*/ + +int PentaElementUnwrap(void* vpthis) +{ + /* though no extended data means no internal pointers to + reestablish, re-null pointers whose values must properly come + from a call to Configure(): */ + + int i; + PentaElement* this = (PentaElement*)vpthis; + + for ( i=0; ipg[i] = NULL; + + this->matice = NULL; + this->matpar = NULL; + + return 1; +} + +/*-------------------------------------------------- + PentaElementMarshall + --------------------------------------------------*/ + +void PentaElementMarshall( char* *pbuffer, void* vpthis, int* size) +{ + /* vpthis for polymorphic function compatibility */ + PentaElement* this = (PentaElement*)vpthis; + + char* buffer=*pbuffer; + + int i, j, size_penta=0; + int ig; + + int num_ints=0, num_doubles=0; + + + /* first marshall the embedded Penta object: ... */ + + PentaMarshall( &buffer, &this->penta, &size_penta); + + + /* ... then marshall the PentaElement member data: */ + + for ( i=0; iinternalgrid_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->matice_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matice_enum, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matice, &buffer, INTEGER_FLAG); num_ints++; + + cellmar( &this->matpar_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matpar_enum, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matpar, &buffer, INTEGER_FLAG); num_ints++; + + *size = size_penta + num_ints*sizeof(int) + num_doubles*sizeof(double); + *pbuffer = buffer; + + return; +} + +/*-------------------------------------------------- + PentaElementDemarshall + --------------------------------------------------*/ + +void* PentaElementDemarshall( DataSet* ds, char* *pbuffer, int machine_flag) +{ + int i, j, ig; + Penta* ppenta=NULL; /* for embedded object */ + + char* buffer=*pbuffer; + PentaElement* this=NULL; + + + /* first demarshall the embedded Penta object: ... */ + + ppenta = PentaDemarshall( NULL, &buffer, machine_flag); + /* create the derived object: */ + this = NewPentaElement( ppenta); + /* and free penta now that it's embedded in the derived object: */ + DeletePenta( (void**)&ppenta); + + + /* ... then demarshall the PentaElement member data: */ + + for ( i=0; iinternalgrid_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->matice_index, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->matice_enum, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->matice, &buffer, INTEGER_FLAG, machine_flag); + + celldem( &this->matpar_index, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->matpar_enum, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->matpar, &buffer, INTEGER_FLAG, machine_flag); + + + /*Set pointers to NULL :*/ + for (i=0;ipg[i]=NULL; + } + this->matice=NULL; + this->matpar=NULL; + + if (ds) PentaElementAddToDataSet( ds, this); + *pbuffer=buffer; + + return (void*)this; +} + +/*-------------------------------------------------- + PentaElementGetID + --------------------------------------------------*/ + +int PentaElementGetID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + PentaElement* this = (PentaElement*)vpthis; + + return PentaGetID( &this->penta); +} + + +/*-------------------------------------------------- + PentaElementSetVirtualFunctions + --------------------------------------------------*/ + +int PentaElementSetVirtualFunctions( void* vp) +{ + + /* VirtualFunctions implemented as a "friend" class: */ + + /* void* for polymorphic function compatibility */ + VirtualFunctions* vf = (VirtualFunctions*)vp; + + if ( VirtualFunctionsInit(vf) ) { + + /* set PentaElement-specific virtual function pointers: */ + + /* general: */ + + vf->Configure = &PentaElementConfigure; + vf->GetName = &PentaElementGetName; + + /* element matrix generation: */ + vf->CreateKMatrix = &PentaElementCreateKMatrix; + + /* loads vector generation: */ + vf->CreatePVector = &PentaElementCreatePVector; + + /* velocity difference generation: */ + vf->CreateDuVector = &PentaElementCreateDuVector; + vf->CreateGradjVectors = &PentaElementCreateGradjVectors; + + /* velocity misfit: */ + vf->Misfit=&PentaElementMisfit; + + return 1; + } + else + return 0; +} + + +/*-------------------------------------------------- + PentaElementAddToDataSet + --------------------------------------------------*/ + +int PentaElementAddToDataSet( DataSet* dataset, PentaElement* 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, PentaElementEnum(), (void*)this, + &DeletePentaElement, &PentaElementEcho, &PentaElementSizeOfInDoubles, + &PentaElementShrinkWrap, &PentaElementUnwrap, + &PentaElementMarshall, + &PentaElementGetID, &PentaElementSetVirtualFunctions); + + return 1; +} + + + +/*-------------------------------------------------- + PentaElementGetName + --------------------------------------------------*/ + +char* PentaElementGetName( void* vpthis) +{ + PentaElement* this = (PentaElement*)vpthis; + + #ifdef _DEBUG_ + _printf_("virtual function: PentaElementGetName\n"); + #endif + + return PentaGetName( &this->penta); +} + +/*-------------------------------------------------- + PentaElementGetNumberOfGrids + --------------------------------------------------*/ + +int PentaElementGetNumberOfGrids( void* vpthis ) +{ + return numgrids; +} + +/*-------------------------------------------------- + PentaElementConfigure + --------------------------------------------------*/ + +int PentaElementConfigure( void* vpthis, int num_datasets, ...) +{ + + /* set up data connections to all necessary, related objects: */ + + PentaElement* this = (PentaElement*)vpthis; + + int i,j, found, foundgrid, noerr=1, first=1; + + /* expected DataSet input: */ + + DataSet* pgpdt=NULL; + DataSet* pgpdtb=NULL; + DataSet* pmpt=NULL; + + /* variable argument list handling: */ + + va_list args; + DataSet* *arglist=NULL; + + #ifdef _DEBUG_ + _printf_("virtual function: PentaElementConfigure\n"); + #endif + + if (this) { + + if (!num_datasets) { + + _printf_(" %s error:\n", __func__); + _printf_(" null input for number of datasets;\n"); + _printf_(" could not establish grid, property and material data 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;ipg[i]=NULL; + } + + /* grid connections based on DataSet of type "gpdt": */ + + #ifndef _PARALLEL_ + found=0; + for ( i=0; ipenta); + + int i; + + for ( i=0; iinternalgrid_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) { + 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, PentaGetID(&this->penta), PentaGetName(&this->penta)); + } + } + } + } + } + else { + + if (first) { + _printf_(" %s error:\n", __func__); + first=0; + } + _printf_(" could not find DataSet of type \"gpdt\", grid object\n"); + _printf_(" pointers not established for penta element id = %d\n", + PentaElementGetID(this)); + noerr = 0; + } + #else //ifdef _PARALLEL_ + + /*Configure the element when we are running on a parallel cluster: in this case, + we have partitioned the elements and the grids correspondingly. Every node has + its own est, 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 + element, 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; ipenta); + + for ( i=0; iinternalgrid_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_indices[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, PentaGetID(&this->penta), PentaGetName(&this->penta)); + } + } + } + } + } + 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 penta element id = %d\n", + PentaElementGetID(this)); + noerr = 0; + } + + #endif + + + /* matice property connections based on DataSet of type "mpt": */ + + found=0; + for ( i=0; ipenta); + + /* as for grids, check to see if a nonzero matice index still + points to the matice entry of interest in mpt, otherwise do + a seek for it: */ + + if ( this->matice_enum && this->matice_index ) { + + + if ( this->matice_enum == DataSetGetEnum( pmpt, this->matice_index) + && + current_matice_id == DataSetGetObjectID( pmpt, this->matice_index) ) { + + /* index is still valid; just get matice entry pointer directly: */ + + this->matice = DataSetGetObjectPtr( pmpt, this->matice_index); + } + } + + /* if we haven't yet "recovered" the matice object pointer, go seek it out: */ + + if ( NULL == this->matice) { + + noerr *= DataSetSeek( &this->matice, &this->matice_index, + pmpt, current_matice_id, 1, MaticeEnum() ); /* 0 to generalise later to Pcomp, etc. ... */ + + if (!noerr) { + if (first) { + _printf_(" %s error:\n", __func__); + first = 0; + } + _printf_(" could not find matice entry with id of %d for element %d of type \"%s\"\n", + current_matice_id, PentaGetID(&this->penta), PentaGetName(&this->penta)); + } + else { + + /* though the answer for now is just MaticeEnum(), for generality go ahead and + ask the question, "what type of matice has been referenced?" */ + + this->matice_enum = DataSetGetEnum( *(arglist+i), this->matice_index); + } + } + + + + int current_matpar_id = PentaGetMaterialParID(&this->penta); + + /* as for grids, check to see if a nonzero matpar index still + points to the matpar entry of interest in mpt, otherwise do + a seek for it: */ + + if ( this->matpar_enum && this->matpar_index ) { + + + if ( this->matpar_enum == DataSetGetEnum( pmpt, this->matpar_index) + && + current_matpar_id == DataSetGetObjectID( pmpt, this->matpar_index) ) { + + /* index is still valid; just get matpar entry pointer directly: */ + + this->matpar = DataSetGetObjectPtr( pmpt, this->matpar_index); + } + } + + /* if we haven't yet "recovered" the matpar object pointer, go seek it out: */ + + if ( NULL == this->matpar) { + + noerr *= DataSetSeek( &this->matpar, &this->matpar_index, + pmpt, current_matpar_id, 1, MatparEnum() ); /* 0 to generalise later to Pcomp, etc. ... */ + + if (!noerr) { + if (first) { + _printf_(" %s error:\n", __func__); + first = 0; + } + _printf_(" could not find matpar entry with id of %d for element %d of type \"%s\"\n", + current_matpar_id, PentaGetID(&this->penta), PentaGetName(&this->penta)); + } + else { + + /* though the answer for now is just MatparEnum(), for generality go ahead and + ask the question, "what type of matpar has been referenced?" */ + + this->matpar_enum = DataSetGetEnum( *(arglist+i), this->matpar_index); + } + } + } + else { + + if (first) { + _printf_(" %s error:\n", __func__); + first=0; + } + _printf_(" could not find DataSet of type \"mpt\"; matice and matpar relationships\n"); + _printf_(" not established for penta element id = %d\n", + PentaElementGetID(this)); + noerr = 0; + } + } + + return noerr; +} +/*-------------------------------------------------- + PentaElementGetGridIDPtr + --------------------------------------------------*/ + +int* PentaElementGetGridIDPtr( void* vpthis) +{ + PentaElement* this = (PentaElement*)vpthis; + + #ifdef _DEBUG_ + _printf_("virtual function: PentaElementGetGridIDPtr\n"); + #endif + + return PentaGetGridIDPtr( &this->penta); +} + + +/*-------------------------------------------------- + PentaElementCreateKMatrix + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrix" +int PentaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){ + + + int noerr=1; + + + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ + noerr=PentaElementCreateKMatrixDiagnosticHoriz( pKe_gg, vpthis, inputs); + } + else if(analysis_type==DiagnosticBaseVertAnalysisEnum()){ + noerr=PentaElementCreateKMatrixDiagnosticBaseVert( pKe_gg, vpthis, inputs); + } + else if(analysis_type==DiagnosticVertAnalysisEnum()){ + noerr=PentaElementCreateKMatrixDiagnosticVert( pKe_gg, vpthis, inputs); + } + else if(analysis_type==ThermalSteadyAnalysisEnum()|| (analysis_type==ThermalTransientAnalysisEnum())){ + noerr=PentaElementCreateKMatrixThermal( pKe_gg, vpthis, inputs); + } + else if(analysis_type==MeltingAnalysisEnum()){ + noerr=PentaElementCreateKMatrixMelting( pKe_gg, vpthis, inputs); + } + else if(analysis_type==PrognosticAnalysisEnum()){ + noerr=PentaElementCreateKMatrixPrognostic( pKe_gg, vpthis, inputs); + } + else{ + _printf_("%s%s%i%s\n",__FUNCT__," error message: analysis type",analysis_type," not supported yet"); + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + PentaElementCreateKMatrixDiagnosticHoriz + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrixDiagnosticHoriz" +int PentaElementCreateKMatrixDiagnosticHoriz( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + + int i,j,il,counter; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + Matice * matice = NULL; + Matpar * matpar = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + /* 3d gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* fourth_gauss_vert_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_l1l2l3l4[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + double gauss_weight; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + + /*drag: */ + double pcoeff,qcoeff; + double rcoeff,scoeff; + double velocity_x,velocity_y,velocity_mag; + + /* material data: */ + double viscosity; //viscosity + double B_param; + double gravity=9.8; + double rho_ice,rho_water; + + /* strain rate: */ + double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + + /* matrices: */ + double B[5][NDOF2*numgrids]; + double Bprime[5][NDOF2*numgrids]; + double L[2][NDOF2*numgrids]; + double D[5][5]={{ 0,0,0,0,0 },{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}}; // material matrix, simple scalar matrix. + double D_scalar; + + double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + double DL_scalar; + + double Ke_gg_gaussian[NDOF2*numgrids][NDOF2*numgrids]; //stiffness matrix evaluated at the gaussian point. + + double Ke_gg_drag_gaussian[NDOF2*3][NDOF2*3]; //stiffness matrix contribution from drag, in 2d + double Jdet; + + /*slope: */ + double slope[NDOF2]; + double slope_magnitude; + + /*input parameters for structural analysis (diagnostic): */ + double* thickness_param=NULL; + double thickness_list[numgrids]; + double thickness; + double* surface_param=NULL; + double surface_list[numgrids]; + double* bed_param=NULL; + double bed_list[numgrids]; + double* flow_law_param=NULL; + double B_list[numgrids]; + double* basal_drag=NULL; + double K_list[numgrids]; + double* velocity=NULL; + double vxvy_list[numgrids][2]; + double* temperature_average_param=NULL; + double temperature_average_list[numgrids]; + double temperature_average; + + //tria element parameters + double alpha2_list[3]; + double alpha2; + + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + + /*Collapsed formulation: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /* The number of g-set dof for this element is 2*6 (vx and vy, horizontal velocities for ice, and 6 grids + * per element: */ + numdof=numgrids*NDOF2; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + matice=this->matice; + matpar=this->matpar; + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the stiffness matrix. */ + + if ((this->penta.collapse==1) && (this->penta.onbed==0)){ + /*This element should be collapsed, but this element is not on the bedrock, therefore all its + * dofs have already been frozen! Do nothing: */ + *pKe_gg=NULL; + return noerr; + } + else if ((this->penta.collapse==1) && (this->penta.onbed==1)){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + *and use its CreateKMatrix functionality to return an elementary stiffness matrix: */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[i]; + tria->s[i]=this->penta.s[i]; + tria->b[i]=this->penta.b[i]; + tria->k[i]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[i]; + triaelement->pg[i]=this->pg[i]; + } + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && PentaElementGetID(this)==ELID){ printf("El id %i Rank %i TriaElement echo in collapsed Penta: \n",ELID,RANK); TriaElementEcho(triaelement); } + #endif + + /*Ok, now triaelement is correctly configured, call on its method CreateKMatrix: */ + noerr=TriaElementCreateKMatrixDiagnosticHoriz( pKe_gg, (void*)triaelement, inputs); + + /*Now delete tria and triaelement: */ + DeleteTriaElement((void**)&triaelement); + DeleteTria((void**)&tria); + + return noerr; + } + else{ + + /*Implement standard penta element: */ + + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /*Recover material: */ + gravity=matpar->g; + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + + /*recover extra inputs from users, at current convergence iteration: */ + basal_drag=ParameterInputsRecover(inputs,"drag"); + thickness_param=ParameterInputsRecover(inputs,"thickness"); + surface_param=ParameterInputsRecover(inputs,"surface"); + bed_param=ParameterInputsRecover(inputs,"bed"); + temperature_average_param=ParameterInputsRecover(inputs,"temperature_average"); + flow_law_param=ParameterInputsRecover(inputs,"B"); + velocity=ParameterInputsRecover(inputs,"velocity"); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=dof; + if(velocity){ + vxvy_list[i][j]=*(velocity+dof); + } + else{ + vxvy_list[i][j]=0; + } + } + dof=structural_dof_list[0]; + if(flow_law_param){ + B_list[i] = *(flow_law_param+dof); + } + if(basal_drag){ + K_list[i]=*(basal_drag+dof); + } + else{ + K_list[i]=this->penta.k[i]; + } + if(thickness_param){ + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->penta.h[i]; + } + if(surface_param){ + surface_list[i]= *(surface_param+dof); + } + else{ + surface_list[i]= this->penta.s[i]; + } + if(bed_param){ + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]= this->penta.b[i]; + } + + if(temperature_average_param) temperature_average_list[i]=*(temperature_average_param+dof); + } + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && PentaElementGetID(this)==ELID){ + printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK); + printf(" rho_ice: %g\n",rho_ice); + printf(" gravity: %g\n",gravity); + printf(" rho_water:%g \n",rho_water); + printf(" Velocity: \n"); + for (i=0;imatice,B_param); + } + //Update material if B is provided. + if(flow_law_param){ + PentaElementGetParameterValue(&B_param, &B_list[0],gauss_l1l2l3l4); + #ifdef _DEBUG_ + _printf_("Flow law parameters: %lf\n", B_param); + #endif + MaticeSetFlowLawParameter(this->matice,B_param); + } + + + /*Get viscosity: */ + MaticeGetViscosity3d(&viscosity, this->matice, &epsilon[0]); + #ifdef _DEBUG_ + _printf_("Element id %i Viscosity: \n", PentaElementGetID(this),viscosity); + #endif + + /*Get B and Bprime matrices: */ + PentaElementGetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3l4); + PentaElementGetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3l4); + + /* Get Jacobian determinant: */ + PentaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4); + + /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + D_scalar=viscosity*gauss_weight*Jdet; + for (i=0;i<5;i++){ + D[i][i]=D_scalar; + } + + + /* Do the triple product tB*D*Bprime: */ + noerr=TripleMultiply( &B[0][0],5,numgrids*NDOF2,1, + &D[0][0],5,5,0, + &Bprime[0][0],5,numgrids*NDOF2,0, + &Ke_gg_gaussian[0][0],0); + + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; inrows; i++){ + for (j=0;jnrows;j++){ + *(Ke_gg->terms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + } //for (ig2=0; ig2penta.onbed) && (this->penta.shelf==0)){ + + /*Build a tria element using the 3 grids of the base of the penta. Then use + * the tria functionality to build a friction stiffness matrix on these 3 + * grids: */ + + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[i]; + tria->s[i]=this->penta.s[i]; + tria->b[i]=this->penta.b[i]; + tria->k[i]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[i]; + triaelement->pg[i]=this->pg[i]; + } + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /*Build alpha2_list used by drag stiffness matrix*/ + if (triaelement->tria.friction_type==2){ + + /*Allocate friction object: */ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=gravity; + friction->rho_ice=rho_ice; + friction->rho_water=rho_water; + friction->K=&K_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=triaelement->tria.p; + friction->q=triaelement->tria.q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igterms+Ke_gg->nrows*i+j)+=Ke_gg_drag_gaussian[i][j]; + } + } + } //for (ig=0; igpenta.onbed) && (this->penta.shelf==0)) + } //else if ((this->tria.collapse==1) && (this->tria.onbed==1)) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&fourth_gauss_vert_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build stiffness matrix."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return 1; + } +} + +/*-------------------------------------------------- + PentaElementCreateKMatrixDiagnosticBaseVert + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrixDiagnosticBaseVert" +int PentaElementCreateKMatrixDiagnosticBaseVert( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int numgrids2d=3; + int numdof; + int* structural_dof_list = NULL; + int dof; + + /* grid data: */ + int cid_list[numgrids2d]; + double xyz_list[numgrids2d][3]; + double T_bg_list[numgrids2d][9]; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + double gauss_weight; + int ig; + + /* matrices: */ + double L[1][3]; + double DL_scalar; + + double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point. + double Jdet; + /*Collapsed formulation: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /* The number of g-set dof for this element is 1*3 (vz on bedrock grids) : */ + numdof=numgrids2d*NDOF1; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*We are dealing with a collapsed formulation on the lower triangle of the penta, + only on the bedrock.:*/ + + if (this->penta.onbed==0){ + *pKe_gg=NULL; + return noerr; + } + else{ + /*Build tria element at the base: */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[i]; + tria->s[i]=this->penta.s[i]; + tria->b[i]=this->penta.b[i]; + tria->k[i]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[i]; + triaelement->pg[i]=this->pg[i]; + } + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, triaelement->pg, numgrids2d, 1 ); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + dof=structural_dof_list[2]; //computing vz + *(Ke_gg->row_indices+i)=dof; + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igterms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + } //for (ig=0; igpenta.onbed==0) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build stiffness matrix."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return 1; + } +} + +/*-------------------------------------------------- + PentaElementCreateKMatrixPrognostic + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrixPrognostic" +int PentaElementCreateKMatrixPrognostic( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int numgrids2d=3; + int numdof; + int* structural_dof_list = NULL; + int dof; + + /* grid data: */ + int cid_list[numgrids2d]; + double xyz_list[numgrids2d][3]; + double T_bg_list[numgrids2d][9]; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + double gauss_weight; + int ig; + + /*parameters: */ + double* velocity=NULL; + double vx_list[numgrids2d]; + double vx; + double vy_list[numgrids2d]; + double vy; + double* dt_param=NULL; + double dt; + + /* matrices: */ + double L[1][3]; + double B[2][numgrids2d]; + double DL_scalar; + double D[2][2]; + + double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point. + double Jdet; + + /* The number of g-set dof for this element is 1*3 (vz on bedrock grids) : */ + numdof=numgrids2d*NDOF1; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*We are dealing with a collapsed formulation on the lower triangle of the penta, + only on the bedrock.:*/ + if (this->penta.onbed==0){ + *pKe_gg=NULL; + return noerr; + } + else{ + + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids2d, 1 ); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /*recover extra inputs from users, at current convergence iteration: */ + dt_param=ParameterInputsRecover(inputs,"dt"); + + if(!dt_param){ + _printf_("%s%s\n",__FUNCT__," error message: missing dt parameter!"); + noerr=0; goto cleanup_and_return; + } + dt=*dt_param; + + if (this->penta.artdiff){ + velocity=ParameterInputsRecover(inputs,"velocity_average"); + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + dof=structural_dof_list[0]; //1st dof holds new thickness + *(Ke_gg->row_indices+i)=dof; + + if(this->penta.artdiff){ + dof=structural_dof_list[0]; + vx_list[i]=*(velocity+dof); + dof=structural_dof_list[0]; + vy_list[i]=*(velocity+dof); + } + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igterms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + } //for (ig=0; igpenta.artdiff){ + for (ig=0; igterms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + + } //for (ig=0; igpenta.artdiff) + + } //if (this->penta.onbed==0) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build stiffness matrix."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return 1; + } +} +/*-------------------------------------------------- + PentaElementCreateKMatrixMelting + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrixMelting" +int PentaElementCreateKMatrixMelting( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + int numgrids2d=3; + int numdof; + int* structural_dof_list = NULL; + int dof; + + /* grid data: */ + int cid_list[numgrids2d]; + double xyz_list[numgrids2d][3]; + double T_bg_list[numgrids2d][9]; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + double gauss_weight; + int ig; + + /* matrices: */ + double L[1][3]; + double DL_scalar; + + double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point, for first 3 grids. + double Jdet; + + /* The number of g-set dof for this element is 1*3 (vz on bedrock grids) : */ + numdof=numgrids2d*NDOF1; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*We are dealing with a collapsed formulation on the lower triangle of the penta, + only on the bedrock.:*/ + + if (this->penta.onbed==0){ + *pKe_gg=NULL; + return noerr; + } + else{ + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids2d, 1 ); //only pick up the first 3 + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + dof=structural_dof_list[0]; //melting is found on the 1st dof + *(Ke_gg->row_indices+i)=dof; + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igterms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + } //for (ig=0; igpenta.onbed==0) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build stiffness matrix."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return 1; + } +} + + +/*-------------------------------------------------- + PentaElementCreateKMatrixDiagnosticVert + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrixDiagnosticVert" +int PentaElementCreateKMatrixDiagnosticVert( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + + int noerr=1; + int i,j; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double xyz_list_tria[3][3]; + double T_bg_list[numgrids][9]; + + /* 3d gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* fourth_gauss_vert_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_l1l2l3l4[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + double gauss_weight; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + + /* surface normal: */ + double x4,y4,z4; + double x5,y5,z5; + double x6,y6,z6; + double v46[3]; + double v56[3]; + double normal[3]; + double norm_normal; + double nz; + + /* matrices: */ + double B[1][numgrids]; + double Bprime[1][numgrids]; + double DL_scalar; + double L[3]; + + double Ke_gg_gaussian[numgrids][numgrids]; + double Ke_gg_vert_gaussian[3][3]; + double Jdet; + /*Collapsed formulation: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /* The number of g-set dof for this element is 2*6 (vx and vy, horizontal velocities for ice, and 6 grids + * per element: */ + numdof=numgrids*NDOF1; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*Implement standard penta element: */ + + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + dof=structural_dof_list[2]; //z degree of freedom + *(Ke_gg->row_indices+i)=dof; + } + + /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + order_area_gauss=2; + num_vert_gauss=2; + + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); + #ifdef _DEBUG_ + for (i=0;inrows; i++){ + for (j=0;jnrows;j++){ + *(Ke_gg->terms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + } //for (ig2=0; ig2penta.onsurface==1){ + /*Build a tria element using the 3 grids at the upper base of the penta: */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[3+i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[3+i]; + tria->s[i]=this->penta.s[3+i]; + tria->b[i]=this->penta.b[3+i]; + tria->k[i]=this->penta.k[3+i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[3+i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[3+i]; + triaelement->pg[i]=this->pg[3+i]; + } + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /*Build xyz_list_tria: */ + xyz_list_tria[0][0]=xyz_list[4][0]; + xyz_list_tria[0][1]=xyz_list[4][1]; + xyz_list_tria[0][2]=xyz_list[4][2]; + + xyz_list_tria[1][0]=xyz_list[5][0]; + xyz_list_tria[1][1]=xyz_list[5][1]; + xyz_list_tria[1][2]=xyz_list[5][2]; + + xyz_list_tria[2][0]=xyz_list[6][0]; + xyz_list_tria[2][1]=xyz_list[6][1]; + xyz_list_tria[2][2]=xyz_list[6][2]; + + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /*Build normal vector to the surface:*/ + + x4=xyz_list[0][0]; + y4=xyz_list[0][1]; + z4=xyz_list[0][2]; + + x5=xyz_list[1][0]; + y5=xyz_list[1][1]; + z5=xyz_list[1][2]; + + x6=xyz_list[2][0]; + y6=xyz_list[2][1]; + z6=xyz_list[2][2]; + + v46[0]=x4-x6; + v46[1]=y4-y6; + v46[2]=z4-z6; + + v56[0]=x5-x6; + v56[1]=y5-y6; + v56[2]=z5-z6; + + normal[0]=(y4-y6)*(z5-z6)-(z4-z6)*(y5-y6); + normal[1]=(z4-z6)*(x5-x6)-(x4-x6)*(z5-z6); + normal[2]=(x4-x6)*(y5-y6)-(y4-y6)*(x5-x6); + + norm_normal=sqrt(pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2)); + nz=1.0/norm_normal*normal[2]; + + /* Start looping on the number of gaussian points: */ + for (ig=0; igterms+Ke_gg->nrows*i+j) = *(Ke_gg->terms+Ke_gg->nrows*i+j)+ Ke_gg_vert_gaussian[i-3][j-3]; + } + } + + } //for (ig=0; igpenta.onsurface==1) + + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&fourth_gauss_vert_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build stiffness matrix."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemMatrixEcho(Ke_gg); + #endif + + + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return 1; + } +} + + +/*-------------------------------------------------- + PentaElementCreateKMatrixThermal + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateKMatrixThermal" +int PentaElementCreateKMatrixThermal( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + + int noerr=1; + int i,j; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + Matice * matice = NULL; + Matpar * matpar = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double xyz_list_tria[3][3]; + double T_bg_list[numgrids][9]; + + /* 3d gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* fourth_gauss_vert_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_l1l2l3l4[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + double gauss_weight; + + /* matrices: */ + double D_scalar; + + double B_conduct[3][numgrids]; + double D_conduct[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; + + double B_advect[3][numgrids]; + double Bprime_advect[3][numgrids]; + double D_advect[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; + + double L_transient[numgrids]; + + double Ke_gg_gaussian_conduct[numgrids][numgrids]; + double Ke_gg_gaussian_advect[numgrids][numgrids]; + double Ke_gg_gaussian_transient[numgrids][numgrids]; + + /* input parameters: */ + double* velocity=NULL; + double vx_list[numgrids]; + double vx; + double vy_list[numgrids]; + double vy; + double vz_list[numgrids]; + double vz; + double* dt_param=NULL; + double dt; + + double Jdet; + + /*materials: */ + double rho_ice,heatcapacity,thermalconductivity; + + /* The number of g-set dof for this element is 1*6 (temperature on all grids) + * per element: */ + numdof=numgrids*NDOF1; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + matpar=this->matpar; + + //recover material parameters + rho_ice=matpar->rho_ice; + heatcapacity=matpar->heatcapacity; + thermalconductivity=matpar->thermalconductivity; + + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /*recover extra inputs from users, at current convergence iteration: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + dt_param=ParameterInputsRecover(inputs,"dt"); + + if((!velocity) || (!dt_param)){ + _printf_("%s%s\n",__FUNCT__," error message: missing velocity or dt parameter!"); + noerr=0; goto cleanup_and_return; + } + dt=*dt_param; + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + dof=structural_dof_list[0]; //1st dof for temperatures + *(Ke_gg->row_indices+i)=dof; + + dof=structural_dof_list[0]; + vx_list[i]=*(velocity+dof); + dof=structural_dof_list[1]; + vy_list[i]=*(velocity+dof); + dof=structural_dof_list[2]; + vz_list[i]=*(velocity+dof); + } + + /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore + get tria gaussian points as well as segment gaussian points. For tria gaussian + points, order of integration is 2, because we need to integrate the product tB*D*B' + which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian + points, same deal, which yields 3 gaussian points.*/ + + order_area_gauss=2; + num_vert_gauss=2; + + GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); + #ifdef _DEBUG_ + for (i=0;ipenta.thermal_steadystate==1){ + D_scalar= gauss_weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity)); + } + else{ + D_scalar= dt*gauss_weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity)); + } + for(i=0;i<3;i++){ + D_conduct[i][i]=D_scalar; + } + + /* Do the triple product tB*D*B': */ + noerr=TripleMultiply( &B_conduct[0][0],3,numgrids,1, + &D_conduct[0][0],3,3,0, + &B_conduct[0][0],3,numgrids,0, + &Ke_gg_gaussian_conduct[0][0],0); + + /* deal with advection: */ + + /*Get B and Bprime matrix: */ + + PentaElementGetB_advect(&B_advect[0][0], &xyz_list[0][0], gauss_l1l2l3l4); + PentaElementGetBprime_advect(&Bprime_advect[0][0], &xyz_list[0][0], gauss_l1l2l3l4); + + /*Get velocities: */ + PentaElementGetParameterValue(&vx, &vx_list[0],gauss_l1l2l3l4); + PentaElementGetParameterValue(&vy, &vy_list[0],gauss_l1l2l3l4); + PentaElementGetParameterValue(&vz, &vz_list[0],gauss_l1l2l3l4); + + /*Build D matrix: */ + if (this->penta.thermal_steadystate==1){ + D_scalar= gauss_weight*Jdet; + } + else{ + D_scalar= dt*gauss_weight*Jdet; + } + D_advect[0][0]=D_scalar*vx; + D_advect[1][1]=D_scalar*vy; + D_advect[2][2]=D_scalar*vz; + + /* Do the triple product tB*D*B': */ + noerr=TripleMultiply( &B_advect[0][0],3,numgrids,1, + &D_advect[0][0],3,3,0, + &Bprime_advect[0][0],3,numgrids,0, + &Ke_gg_gaussian_advect[0][0],0); + + /*Deal with transient: */ + if(this->penta.thermal_steadystate==1){ + /*steady state, do nothing: */ + for(i=0;inrows; i++){ + for (j=0;jnrows;j++){ + *(Ke_gg->terms+Ke_gg->nrows*i+j)+= Ke_gg_gaussian_conduct[i][j]+Ke_gg_gaussian_advect[i][j]+Ke_gg_gaussian_transient[i][j]; + } + } + + + } //for (ig2=0; ig2matice; + matpar=this->matpar; + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the load vector. */ + + if ((this->penta.collapse==1) && (this->penta.onbed==0)){ + /*This element should be collapsed, but this element is not on the bedrock, therefore all its + * dofs have already been frozen! Do nothing: */ + *ppe_g=NULL; + return noerr; + } + else if ((this->penta.collapse==1) && (this->penta.onbed==1)){ + + /*This element should be collapsed into a tria element at its base. Create this tria element, + *and use its CreatePVector functionality to return an elementary stiffness matrix: */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[i]; + tria->s[i]=this->penta.s[i]; + tria->b[i]=this->penta.b[i]; + tria->k[i]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[i]; + triaelement->pg[i]=this->pg[i]; + } + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + /*Ok, now triaelement is correctly configured, call on its method CreateKMatrix: */ + noerr=TriaElementCreatePVectorDiagnosticHoriz( ppe_g, (void*)triaelement, inputs); + + /*Now delete tria and triaelement: */ + DeleteTriaElement((void**)&triaelement); + DeleteTria((void**)&tria); + + return noerr; + } + else{ + + /*Implement standard penta element: */ + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids per element: */ + numdof=numgrids*NDOF2; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* recover material parameters: */ + gravity=matpar->g; + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + + /* recover input parameters: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + surface_param=ParameterInputsRecover(inputs,"surface"); + bed_param=ParameterInputsRecover(inputs,"bed"); + melting_param=ParameterInputsRecover(inputs,"melting"); + basal_drag_param=ParameterInputsRecover(inputs,"drag"); + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids, 1 ); + + /* Build the row index vector : */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=structural_dof_list[j]; + } + if(thickness_param){ + dof=structural_dof_list[0]; + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->penta.h[i]; + } + + if(surface_param){ + dof=structural_dof_list[0]; + surface_list[i]=*(surface_param+dof); + } + else{ + surface_list[i]=this->penta.s[i]; + } + if(bed_param){ + dof=structural_dof_list[0]; + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]=this->penta.b[i]; + } + if(melting_param){ + dof=structural_dof_list[0]; + melting_list[i]=*(melting_param+dof); + } + else{ + melting_list[i]=this->penta.b[i]; + } + if(basal_drag_param){ + dof=structural_dof_list[0]; + K_list[i]=*(basal_drag_param+dof); + } + else{ + K_list[i]=this->penta.k[i]; + } + } //for (i=0;ipenta.shelf==0) && (this->penta.friction_type==1)){ + PentaElementGetParameterValue(&plastic_stress, &K_list[0],gauss_l1l2l3l4); + #ifdef _DEBUG_ + _printf_("Plastic stress: %lf\n", plastic_stress); + #endif + } + + /* Get Jacobian determinant: */ + PentaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4); + #ifdef _DEBUG_ + _printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet); + #endif + + /*Get nodal functions: */ + PentaElementGetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4); + + /*Compute driving stress: */ + driving_stress_baseline=rho_ice*gravity; + #ifdef _DEBUG_ + _printf_("Element id %i driving_stress_baseline: %lf\n",PentaElementGetID(this),driving_stress_baseline); + #endif + + /*Build pe_g_gaussian vector: */ + if(this->penta.friction_type==1){ + for (i=0;inrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + } //for (ig2=0; ig2penta.collapse==1) && (this->penta.onbed==1)) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&fourth_gauss_vert_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } +} + +/*-------------------------------------------------- + PentaElementCreatePVectorDiagnosticBaseVert + --------------------------------------------------*/ + +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreatePVectorDiagnosticBaseVert" + +int PentaElementCreatePVectorDiagnosticBaseVert( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + Matice * matice = NULL; + Matpar * matpar = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + int numgrids2d=3; + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* parameters: */ + double rho_ice; + double rho_water; + double di; + double slope[NDOF2]; + double dbdx,dbdy; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + double gauss_weight; + int ig; + + /* grid data: */ + int cid_list[numgrids2d]; + double xyz_list[numgrids2d][3]; + double T_bg_list[numgrids2d][9]; + + /* Jacobian: */ + double Jdet; + + /* matrices: */ + double L[numgrids2d]; + double HU[numgrids2d]; + double dHU[2]; + double HV[numgrids2d]; + double dHV[2]; + double divHVel; + + /*element vector at the gaussian points: */ + double pe_g_gaussian[numgrids2d*NDOF1]; + + /*input parameters for structural analysis (diagnostic): */ + double* thickness_param=NULL; + double thickness_list[numgrids2d]; + double thickness; + double* bed_param=NULL; + double bed_list[numgrids2d]; + double* melting_param=NULL; + double melting_list[numgrids2d]; + double melting; + double* accumulation_param=NULL; + double accumulation_list[numgrids2d]; + double accumulation; + double* velocity_param=NULL; + double vx_list[numgrids2d]; + double vy_list[numgrids2d]; + double vx,vy; + double* velocity_average_param=NULL; + double vx_average_list[numgrids2d]; + double vy_average_list[numgrids2d]; + + + /*Collapsed formulation: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + matice=this->matice; + matpar=this->matpar; + + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the load vector. */ + + if (this->penta.onbed==0){ + *ppe_g=NULL; + return noerr; + } + else{ + + /*Implement standard penta element: */ + + /* The number of g-set dof for this element is 3*1 (vz for 3 base grids: */ + numdof=numgrids2d*NDOF1; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* recover material parameters: */ + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + di=rho_ice/rho_water; + + /* recover input parameters: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + bed_param=ParameterInputsRecover(inputs,"bed"); + melting_param=ParameterInputsRecover(inputs,"melting"); + accumulation_param=ParameterInputsRecover(inputs,"accumulation"); + velocity_param=ParameterInputsRecover(inputs,"velocity"); + velocity_average_param=ParameterInputsRecover(inputs,"velocity_average"); + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids2d, 1 ); //just the first 3 grids + + /* Initialize row vector and parameter lists: */ + for (i=0;ipg[i]); + dof=structural_dof_list[2]; //dof for vz + *(pe_g->row_indices+i)=dof; + + /* Recover vx: */ + dof=structural_dof_list[0]; //dof for vx + if(velocity_param){ + vx_list[i]=*(velocity_param+dof); + } + else{ + vx_list[i]=0; + } + if(velocity_average_param){ + vx_average_list[i]=*(velocity_average_param+dof); + } + else{ + vx_average_list[i]=0; + } + /* Recover vy: */ + dof=structural_dof_list[1]; //dof for vy + if(velocity_param){ + vy_list[i]=*(velocity_param+dof); + } + else{ + vy_list[i]=0; + } + if(velocity_average_param){ + vy_average_list[i]=*(velocity_average_param+dof); + } + else{ + vy_average_list[i]=0; + } + + /*Recover other parameters on 1 dof : */ + if(thickness_param){ + dof=structural_dof_list[0]; //dof for vx, where thickness is + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->penta.h[i]; + } + if(bed_param){ + dof=structural_dof_list[0]; + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]=this->penta.b[i]; + } + if(melting_param){ + dof=structural_dof_list[0]; + melting_list[i]=*(melting_param+dof); + } + else{ + melting_list[i]=this->penta.melting[i]; + } + if(accumulation_param){ + dof=structural_dof_list[0]; + accumulation_list[i]=*(accumulation_param+dof); + } + else{ + accumulation_list[i]=this->penta.accumulation[i]; + } + } //for (i=0;iname,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[i]; + tria->s[i]=this->penta.s[i]; + tria->b[i]=this->penta.b[i]; + tria->k[i]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[i]; + triaelement->pg[i]=this->pg[i]; + } + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /*For icesheets: */ + /* Start looping on the number of gaussian points: */ + for (ig=0; ignrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + } //for (ig=0; ig we add a corrective term + to the icesheet term.:*/ + if (this->penta.shelf){ + + /*build the vector thickness*velocity_average:*/ + for(i=0;inrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + } //for (ig=0; igpenta.shelf) + + /*Delete Tria and TriaElement objects: */ + DeleteTriaElement((void**)&triaelement); + DeleteTria((void**)&tria); + } //if (this->penta.onbed==0) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } +} + +/*-------------------------------------------------- + PentaElementCreatePVectorPrognostic + --------------------------------------------------*/ + +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreatePVectorPrognostic" + +int PentaElementCreatePVectorPrognostic( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + int numgrids2d=3; + int numdof; + int* structural_dof_list = NULL; + int dof; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + double gauss_weight; + int ig; + + /* grid data: */ + int cid_list[numgrids2d]; + double xyz_list[numgrids2d][3]; + double T_bg_list[numgrids2d][9]; + + /* Jacobian: */ + double Jdet; + + /* matrices: */ + + /*element vector at the gaussian points: */ + double pe_g_gaussian[numgrids2d*NDOF1]; + + /*input parameters for structural analysis (diagnostic): */ + double* basal_vertical_velocity_param=NULL; + double wb_list[numgrids2d]; + double wb; + double* surface_vertical_velocity_param=NULL; + double ws_list[numgrids2d]; + double ws; + double* accumulation_param=NULL; + double accumulation_list[numgrids2d]; + double accumulation; + double* melting_param=NULL; + double melting_list[numgrids2d]; + double melting; + double* thickness_param=NULL; + double thickness_list[numgrids2d]; + double thickness; + double* dt_param=NULL; + double dt; + + /*Matrices: */ + double L[numgrids2d]; + + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the load vector. */ + + if (this->penta.onbed==0){ + *ppe_g=NULL; + return noerr; + } + else{ + + /*Implement standard penta element: */ + + /* The number of g-set dof for this element is 3*1 (vz for 3 base grids: */ + numdof=numgrids2d*NDOF1; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* recover input parameters: */ + accumulation_param=ParameterInputsRecover(inputs,"accumulation"); + melting_param=ParameterInputsRecover(inputs,"melting"); + basal_vertical_velocity_param=ParameterInputsRecover(inputs,"basal_vertical_velocity"); + surface_vertical_velocity_param=ParameterInputsRecover(inputs,"surface_vertical_velocity"); + thickness_param=ParameterInputsRecover(inputs,"thickness"); + dt_param=ParameterInputsRecover(inputs,"dt"); + + if(!accumulation_param || !melting_param || !basal_vertical_velocity_param || !surface_vertical_velocity_param || !thickness_param || !dt_param){ + _printf_("%s%s\n",__FUNCT__," error message: missing input parameters!"); + _printf_("Here is the list of pointers: \n"); + _printf_("accumulation: %p\n",accumulation_param); + _printf_("melting: %p\n",melting_param); + _printf_("basal_vertical_velocity: %p\n",basal_vertical_velocity_param); + _printf_("surface_vertical_velocity: %p\n",surface_vertical_velocity_param); + _printf_("thickness: %p\n",thickness_param); + _printf_("dt: %p\n",dt_param); + noerr=0; goto cleanup_and_return; + } + dt=*dt_param; + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids2d, 1 ); //just the first 3 grids + + /* Initialize row vector and parameter lists: */ + for (i=0;ipg[i]); + dof=structural_dof_list[0]; //dof for thickness + *(pe_g->row_indices+i)=dof; + + /* Recover parameters on 1st dof: */ + dof=structural_dof_list[0]; + wb_list[i]=*(basal_vertical_velocity_param+dof); + ws_list[i]=*(surface_vertical_velocity_param+dof); + accumulation_list[i]=*(accumulation_param+dof); + melting_list[i]=*(melting_param+dof); + thickness_list[i]=*(thickness_param+dof); + } //for (i=0;inrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + } //for (ig=0; igpenta.onbed==0) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } +} +/*-------------------------------------------------- + PentaElementCreatePVectorMelting + --------------------------------------------------*/ + +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreatePVectorMelting" + +int PentaElementCreatePVectorMelting( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + Matice * matice = NULL; + Matice * matpar = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + int numgrids2d=3; + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* parameters: */ + double rho_ice; + double rho_water; + double di; + double slope[NDOF2]; + double dbdx,dbdy; + + /* 2d gaussian point: */ + int num_gauss2d; + double* first_gauss_area_coord2d = NULL; + double* second_gauss_area_coord2d = NULL; + double* third_gauss_area_coord2d = NULL; + double* gauss_weights2d=NULL; + double gauss_l1l2l3[3]; + double gauss_weight; + int ig; + + /* grid data: */ + int cid_list[numgrids2d]; + double xyz_list[numgrids2d][3]; + double T_bg_list[numgrids2d][9]; + + /* Jacobian: */ + double Jdet; + + /* matrices: */ + double L[numgrids2d]; + double HU[numgrids2d]; + double dHU[2]; + double HV[numgrids2d]; + double dHV[2]; + double divHVel; + + /*element vector at the gaussian points: */ + double pe_g_gaussian[numgrids2d*NDOF1]; + + /*input parameters for structural analysis (diagnostic): */ + double* beta_param=NULL; + double beta; + double* latentheat_param=NULL; + double latentheat; + double* heatcapacity_param=NULL; + double heatcapacity; + double* meltingpoint_param=NULL; + double meltingpoint; + double* pressure_param=NULL; + double pressure_list[numgrids2d]; + double* temperature_param=NULL; + double temperature_list[numgrids2d]; + double* dt_param=NULL; + double dt; + double* penalty_param=NULL; + double penalty; + + double Tpmp_list[numgrids2d]; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the + bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build + the load vector. */ + + if (this->penta.onbed==0){ + *ppe_g=NULL; + return noerr; + } + else{ + + /*Implement standard penta element: */ + + /* The number of g-set dof for this element is 3*1 (vz for 3 base grids: */ + numdof=numgrids2d*NDOF1; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* recover input parameters: */ + beta_param=ParameterInputsRecover(inputs,"beta"); + latentheat_param=ParameterInputsRecover(inputs,"latentheat"); + heatcapacity_param=ParameterInputsRecover(inputs,"heatcapacity"); + meltingpoint_param=ParameterInputsRecover(inputs,"meltingpoint"); + pressure_param=ParameterInputsRecover(inputs,"pressure"); + temperature_param=ParameterInputsRecover(inputs,"temperature"); + dt_param=ParameterInputsRecover(inputs,"dt"); + penalty_param=ParameterInputsRecover(inputs,"penalty"); + + if(!beta_param || !latentheat_param || !heatcapacity_param || !meltingpoint_param || !pressure_param || !temperature_param || !dt_param || !penalty_param){ + _printf_("%s%s\n",__FUNCT__," error message: missing input parameters!"); + _printf_("Here is the list of pointers: \n"); + _printf_("beta: %p\n",beta_param); + _printf_("latentheat: %p\n",latentheat_param); + _printf_("heatcapacity: %p\n",heatcapacity_param); + _printf_("meltingpoint: %p\n",meltingpoint_param); + _printf_("pressure: %p\n",pressure_param); + _printf_("temperature: %p\n",temperature_param); + _printf_("dt: %p\n",dt_param); + _printf_("penalty: %p\n",penalty_param); + noerr=0; goto cleanup_and_return; + } + beta=*beta_param; + latentheat=*latentheat_param; + heatcapacity=*heatcapacity_param; + meltingpoint=*meltingpoint_param; + dt=*dt_param; + penalty=*penalty_param; + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids2d, 1 ); //just the first 3 grids + + /* Initialize row vector and parameter lists: */ + for (i=0;ipg[i]); + dof=structural_dof_list[0]; //melting + *(pe_g->row_indices+i)=dof; + + /*Recover other parameters on 1 dof : */ + temperature_list[i]=*(temperature_param+dof); + pressure_list[i]=*(pressure_param+dof); + } //for (i=0;inrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + } //if (this->penta.onbed==0) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } +} +/*-------------------------------------------------- + PentaElementCreatePVectorDiagnosticVert + --------------------------------------------------*/ + +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreatePVectorDiagnosticVert" + +int PentaElementCreatePVectorDiagnosticVert( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* fourth_gauss_vert_coord = NULL; + double* area_gauss_weights = NULL; + double* vert_gauss_weights = NULL; + double gauss_l1l2l3l4[4]; + int order_area_gauss; + int num_vert_gauss; + int num_area_gauss; + int ig1,ig2; + double gauss_weight1,gauss_weight2; + double gauss_weight; + + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3l4l5l6[6]; + + /*element vector at the gaussian points: */ + double pe_g_gaussian[numgrids]; + + /*input parameters for structural analysis (diagnostic): */ + double* velocity=NULL; + double vx_list[numgrids]; + double vy_list[numgrids]; + double du[3]; + double dv[3]; + double dudx,dvdy; + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids per element: */ + numdof=numgrids*NDOF1; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* recover input parameters: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids, 1 ); + + /* Build the row index vector : */ + for (i=0;ipg[i]); + + //recover dofs for the load vector + dof=structural_dof_list[2]; //vz dof + *(pe_g->row_indices+i)=dof; + + //recover vx + dof=structural_dof_list[0]; //vx dof + if(velocity){ + vx_list[i]=*(velocity+dof); + } + else{ + vx_list[i]=0; + } + //recover vy + dof=structural_dof_list[1]; //vy dof + if(velocity){ + vy_list[i]=*(velocity+dof); + } + else{ + vy_list[i]=0; + } + } //for (i=0;inrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + + } //for (ig2=0; ig2matice; + matpar=this->matpar; + + /*Recover material: */ + gravity=matpar->g; + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + heatcapacity=matpar->heatcapacity; + + /* The number of g-set dof for this element is 6*1 (temperature on all grids):*/ + numdof=numgrids*NDOF1; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids, 1 ); + + /* recover input parameters: */ + velocity=ParameterInputsRecover(inputs,"velocity"); + flow_law_param=ParameterInputsRecover(inputs,"B"); + basal_drag=ParameterInputsRecover(inputs,"drag"); + thickness_param=ParameterInputsRecover(inputs,"thickness"); + temperature_param=ParameterInputsRecover(inputs,"temperature"); + surface_param=ParameterInputsRecover(inputs,"surface"); + bed_param=ParameterInputsRecover(inputs,"bed"); + dt_param=ParameterInputsRecover(inputs,"dt"); + geothermalflux_param=ParameterInputsRecover(inputs,"geothermalflux"); + + if((!velocity) || (!dt_param) ){ + _printf_("%s%s\n",__FUNCT__," error message: missing velocity or dt or temperature parameter!"); + noerr=0; goto cleanup_and_return; + } + dt=*dt_param; + + if(this->penta.thermal_steadystate==0){ + if(!temperature_param){ + _printf_("%s%s\n",__FUNCT__," error message: missing temperature for trnasient analysis!"); + } + } + + /* Build the row index vector : */ + for (i=0;ipg[i]); + + //recover dofs for the load vector + dof=structural_dof_list[0]; //temperature dof + *(pe_g->row_indices+i)=dof; + + //recover vxvyvz_list: */ + for(j=0;jpenta.k[i]; + } + if(temperature_param){ + temperature_list[i]=*(temperature_param+dof); + } + if(thickness_param){ + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->penta.h[i]; + } + if(surface_param){ + surface_list[i]= *(surface_param+dof); + } + else{ + surface_list[i]= this->penta.s[i]; + } + if(bed_param){ + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]= this->penta.b[i]; + } + if(geothermalflux_param){ + geothermalflux_list[i]=*(geothermalflux_param+dof); + } + else{ + geothermalflux_list[i]= this->penta.geothermalflux[i]; + } + } //for (i=0;imatice,B_param); + } + + /*Build deformational heating: */ + PentaElementGetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4); + MaticeGetViscosity3d(&viscosity, this->matice, &epsilon[0]); + + //Build the strain rate matrix: */ + epsilon_matrix[0]=epsilon[0]; + epsilon_matrix[1]=epsilon[1]; + epsilon_matrix[2]=epsilon[2]; + epsilon_matrix[3]=epsilon[3]; + epsilon_matrix[4]=epsilon[4]; + epsilon_matrix[5]=-epsilon[0]-epsilon[1]; //incompressibility! + + epsilon_e=effective_value(epsilon_matrix,6,1); + + /*compute source term: */ + phi=2*viscosity*pow(epsilon_e,2); + + /* Get Jacobian determinant: */ + PentaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4); + #ifdef _DEBUG_ + _printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet); + #endif + + /*Get nodal functions: */ + PentaElementGetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4); + + /*Build deformation pe_g_gaussian vector: */ + if(this->penta.thermal_steadystate==1){ + for (i=0;ipenta.thermal_steadystate==1){ + for (i=0;inrows; i++){ + *(pe_g->terms+i)+= *(pe_g_gaussian_deformational+i) + *(pe_g_gaussian_transient+i); + } + + } //for (ig2=0; ig2penta.onbed) && (this->penta.shelf==0)){ + + //Get the basal friction tau_b.u_b=alpha2*ub^2 + if (this->penta.friction_type==2){ + + /*Allocate friction object: */ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=gravity; + friction->rho_ice=rho_ice; + friction->rho_water=rho_water; + friction->K=&K_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=this->penta.p; + friction->q=this->penta.q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + + /*Get velocity magnitude and compute basal friction: */ + for(i=0;i<3;i++){ + vel_base_squared[i]=pow(vxvy_list[i][0],2)+pow(vxvy_list[i][1],2); + basal_friction_list[i]=alpha2_list[i]*vel_base_squared[i]; + } + } //if (this->penta.friction_type==2) + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss2d, &first_gauss_area_coord2d, &second_gauss_area_coord2d, &third_gauss_area_coord2d, &gauss_weights2d, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igpenta.thermal_steadystate){ + for(i=0;i<3;i++){ + pe_geo_gaussian[i]=gauss_weight*Jdet*(basal_friction+G_g)/(heatcapacity*rho_ice)*l1l2l3[i]; + } + } + else{ + for(i=0;i<3;i++){ + pe_geo_gaussian[i]=dt*gauss_weight*Jdet*(basal_friction+G_g)/(heatcapacity*rho_ice)*l1l2l3[i]; + } + } //if(this->penta.thermal_steadystate) + + /*Add gauss load to global load vector: */ + for( i=0; i<3; i++){ //only the first 3 grids. + *(pe_g->terms+i)+= pe_geo_gaussian[i]; + } + + } //for (ig=0; igpenta.onbed) && (this->penta.shelf==0)) + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&fourth_gauss_vert_coord); + xfree((void**)&area_gauss_weights); + xfree((void**)&vert_gauss_weights); + xfree((void**)&first_gauss_area_coord2d); + xfree((void**)&second_gauss_area_coord2d); + xfree((void**)&third_gauss_area_coord2d); + xfree((void**)&gauss_weights2d); + + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } +} +/*-------------------------------------------------- + PentaElementGetParameterDerivativeValue + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementGetParameterDerivativeValue" +int PentaElementGetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){ + + /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_l1l2l3l4: + * dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx; + * dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy; + * dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz; + * + * p is a vector of size 3x1 already allocated. + */ + + int noerr=1; + double dh1dh2dh3dh4dh5dh6_basic[NDOF3][numgrids]; + + /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */ + noerr=PentaElementGetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4); + if(!noerr)goto cleanup_and_return; + + *(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5]; +; + *(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5]; + + *(p+2)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[2][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[2][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[2][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[2][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[2][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[2][5]; + + cleanup_and_return: + if(!noerr){ + _printf_("%s%s\n",__FUNCT__," error message: could not return parameter derivative vector."); + return 0; + } + else{ + return 1; + } +} + +/*-------------------------------------------------- + PentaElementCreateDuVector + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateDuVector" +int PentaElementCreateDuVector( ElemVector* *pdue_g, void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ + + int noerr=1; + int i; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /*surface tria element: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*Bail out if this element does not touch the surface: */ + if (!this->penta.onsurface){ + *pdue_g=NULL; + } + else{ + + /*use the TriaElementCreateDuVector routine, on a TriaElement which is made of the 3 upper grids at the surface of the ice sheet. : */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=3;i<6;i++){ + tria->g[i-3]=this->penta.g[i]; //we take the last 3 grids of the penta to build the tria. + tria->h[i-3]=this->penta.h[i]; + tria->s[i-3]=this->penta.s[i]; + tria->b[i-3]=this->penta.b[i]; + tria->k[i-3]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=3;i<6;i++){ + triaelement->internalgrid_indices[i-3]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i-3]=this->internalgrid_indicesb[i]; + triaelement->pg[i-3]=this->pg[i]; + } + + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /*Ok, now triaelement is correctly configured, call on its method CreateDuVector: */ + noerr=TriaElementCreateDuVector( pdue_g, (void*)triaelement, velocity, obs_velocity, inputs, analysis_type); + + /*Now delete tria and triaelement: */ + DeleteTriaElement((void**)&triaelement); + DeleteTria((void**)&tria); + + } + return noerr; +} + +/*-------------------------------------------------- + PentaElementCreateGradjVectors + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateGradVectors" +int PentaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){ + + int noerr=1; + + if strcasecmp_eq(control_type,"drag"){ + noerr=PentaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + } + else if strcasecmp_eq(control_type,"B"){ + noerr=PentaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + } + else{ + printf("%s%s\n",__FUNCT__," error message: control type not yet supported!"); + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + PentaElementMisfit + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementMisfit" +int PentaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ + + int noerr=1; + int i; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /*surface tria element: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*Bail out if this element does not touch the surface: */ + if (!this->penta.onsurface){ + *pJelem=0; + } + else{ + + /*use the TriaElementCreateDuVector routine, on a TriaElement which is made of the 3 upper grids at the surface of the ice sheet. : */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=3;i<6;i++){ + tria->g[i-3]=this->penta.g[i]; //we take the last 3 grids of the penta to build the tria. + tria->h[i-3]=this->penta.h[i]; + tria->s[i-3]=this->penta.s[i]; + tria->b[i-3]=this->penta.b[i]; + tria->k[i-3]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=3;i<6;i++){ + triaelement->internalgrid_indices[i-3]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i-3]=this->internalgrid_indicesb[i]; + triaelement->pg[i-3]=this->pg[i]; + } + + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /*Ok, now triaelement is correctly configured, call on its method Misfit: */ + noerr=TriaElementMisfit(pJelem,(void*)triaelement, velocity, obs_velocity, inputs, analysis_type); + + /*Now delete tria and triaelement: */ + DeleteTriaElement((void**)&triaelement); + DeleteTria((void**)&tria); + + } + return noerr; +} + +/*-------------------------------------------------- + PentaElementCreateGradjVectorsDrag + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateGradVectorsDrag" +int PentaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ + + int noerr=1; + int i; + + /* vpthis for polymorphic function compatibility */ + PentaElement* this = NULL; + + /*surface tria element: */ + Tria* tria=NULL; + TriaElement* triaelement=NULL; + + /*Some pointer intialization: */ + this = (PentaElement*)vpthis; + + /*Bail out if this element does not touch the bedrock: */ + if (!this->penta.onbed){ + *pgradje_g=NULL; + } + else{ + + /*use the TriaElementCreateDuVector routine, on a TriaElement which is made of the 3 lower grids at the bedrock interface with the ice sheet: */ + tria=NewTria(); + TriaInit(tria); + + strcpy(tria->name,"CTRICE3"); + tria->eid=this->penta.eid; + tria->mid=this->penta.mid; + for(i=0;i<3;i++){ + tria->g[i]=this->penta.g[i]; //we take the first 3 grids of the penta to build the tria. + tria->h[i]=this->penta.h[i]; + tria->s[i]=this->penta.s[i]; + tria->b[i]=this->penta.b[i]; + tria->k[i]=this->penta.k[i]; + } + + /*diverse:*/ + tria->theta_type=0; + tria->theta.angle=0; + tria->zoffs=0; + + //friction + tria->p=this->penta.p; + tria->q=this->penta.q; + tria->shelf=this->penta.shelf; + tria->friction_type=this->penta.friction_type; + + /*type of fitting? : */ + tria->fit=this->penta.fit; + tria->meanvel=this->penta.meanvel; + tria->epsvel=this->penta.epsvel; + + /*diverse:*/ + tria->acceleration=0; + + /*Now that the tria object is created, spawn an element out of it: */ + triaelement = NewTriaElement( tria); + + /*Transfer configuration of penta to configuration of tria: */ + for(i=0;i<3;i++){ + triaelement->internalgrid_indices[i]=this->internalgrid_indices[i]; + triaelement->internalgrid_indicesb[i]=this->internalgrid_indicesb[i]; + triaelement->pg[i]=this->pg[i]; + } + + triaelement->matice_index=this->matice_index; + triaelement->matice_enum=this->matice_enum; + triaelement->matice=this->matice; + + triaelement->matpar_index=this->matpar_index; + triaelement->matpar_enum=this->matpar_enum; + triaelement->matpar=this->matpar; + + + /*Ok, now triaelement is correctly configured, call on its method CreateGradjVectorsdrag: */ + noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, (void*)triaelement,velocity, adjoint, inputs,analysis_type); + + /*Now delete tria and triaelement: */ + DeleteTriaElement((void**)&triaelement); + DeleteTria((void**)&tria); + + } + return noerr; +} + + +/*-------------------------------------------------- + PentaElementCreateGradjVectorsB + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "PentaElementCreateGradjVectorsB" +int PentaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ + + printf("Not supported yet!\n"); + return 0; +} + + + + +#undef NDOF1 +#undef NDOF2 +#undef NDOF3 +#undef numgrids Index: /issm/trunk/src/c/objects/cielo/TriaElement.c =================================================================== --- /issm/trunk/src/c/objects/cielo/TriaElement.c (revision 216) +++ /issm/trunk/src/c/objects/cielo/TriaElement.c (revision 216) @@ -0,0 +1,2947 @@ +/* + TriaElement.c +*/ + +#include "../../include/cielo.h" + +#define NDOF1 1 +#define NDOF2 2 +#define numgrids 3 + +/*For debugging purposes: */ +#define ELID 141 //element for which to print debug statements +#define RANK 2 //rank of cpu for which to print debug statements. +//#define _DEBUGELEMENTS_ +//#define _DEBUGGAUSS_ + +/*-------------------------------------------------- + NewTriaElement + --------------------------------------------------*/ + +TriaElement* NewTriaElement( void* vptria) +{ + + Tria* tria = (Tria*)vptria; + TriaElement* this = NULL; + int i,j; + + this=(TriaElement*)xcalloc(1,sizeof(TriaElement)); /* guaranteed NULL, 0 initialisation */ + + /* include the "base class" in the "derived" object: */ + memcpy( &this->tria, tria, sizeof(Tria)); + + /* since zero could be a valid index, set all DataSet index + "pointers" to UNDEF: */ + + for ( i=0; i<3; i++){ + this->internalgrid_indices[i] = UNDEF; + this->internalgrid_indicesb[i] = UNDEF; + this->pg[i] = NULL; + } + + this->property_index=UNDEF; + this->property_enum=UNDEF; + this->property=NULL; + + + this->matice_index=UNDEF; + this->matice_enum=UNDEF; + this->matice=NULL; + + this->matpar_index=UNDEF; + this->matpar_enum=UNDEF; + this->matpar=NULL; + + + return this; +} + + +/*-------------------------------------------------- + DeleteTriaElement + --------------------------------------------------*/ + +void DeleteTriaElement(void* *vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + TriaElement* *pthis = (TriaElement**)vpthis; + TriaElement* this = *pthis; + + #ifdef _DEBUG_ + printf(" DeleteTriaElement message:\n"); + printf(" deleting TriaElement with eid = %d\n", this->tria.eid); + #endif + xfree(vpthis); +} + + +/*-------------------------------------------------- + TriaElementEcho + --------------------------------------------------*/ + +void TriaElementEcho(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + TriaElement* this = (TriaElement*)vpthis; + + int i,j; + + printf("\nTriaElement echo:\n"); + printf("-------------------\n"); + printf(" base Tria object:\n"); + TriaEcho( &this->tria); + printf("\n"); + + printf(" element points to following grids (pointer-based ops):\n"); + + printf(" local grid number corresponding grid id\n"); + printf(" ----------------- ---------------------\n"); + + for ( i=0; i<3; i++) { + if ( this->pg[i]) + printf(" %8d %8d\n", i, InternalGridGetID( this->pg[i])); + } + + printf(" element references the following:\n"); + + if ( this->property) printf(" property entry of enum type %d\n", this->property_enum); + else printf(" property reference is NULL\n"); + if ( this->matice) printf(" matice entry of enum type %d\n", this->matice_enum); + else printf(" matice reference is NULL\n"); + if ( this->matpar) printf(" matpar entry of enum type %d\n", this->matpar_enum); + else printf(" matpar reference is NULL\n"); +} + + +/*-------------------------------------------------- + TriaElementSizeOfInDoubles + --------------------------------------------------*/ + +int TriaElementSizeOfInDoubles(void* vpthis) +{ + /* vpthis for polymorphic function compatibility */ + + TriaElement* this = (TriaElement*)vpthis; + int struct_doubles=0; + + struct_doubles = + (int)ceil((float)sizeof(*this)/(float)sizeof(double)); + + #ifdef _DEBUG_ + printf("\n TriaElementSizeOfInDoubles diagnostics:\n"); + printf(" returning %d doubles\n", struct_doubles); + #endif + + return struct_doubles; +} + +/*-------------------------------------------------- + TriaElementShrinkWrap + --------------------------------------------------*/ + +void TriaElementShrinkWrap(void* to, void* vpthis) +{ + /* no "extended" data means a simple copy will do: */ + + memcpy( to, vpthis, TriaElementSizeOfInDoubles(vpthis)*sizeof(double)); +} + + +/*-------------------------------------------------- + TriaElementUnwrap + --------------------------------------------------*/ + +int TriaElementUnwrap(void* vpthis) +{ + /* though no extended data means no internal pointers to + reestablish, re-null pointers whose values must properly come + from a call to Configure(): */ + + TriaElement* this = (TriaElement*)vpthis; + int i; + + for ( i=0; i<3; i++) this->pg[i] = NULL; + this->matice = NULL; + this->matpar = NULL; + this->property = NULL; + return 1; +} + +/*-------------------------------------------------- + TriaElementMarshall + --------------------------------------------------*/ + +void TriaElementMarshall( char* *pbuffer, void* vpthis, int* size) +{ + /* vpthis for polymorphic function compatibility */ + TriaElement* this = (TriaElement*)vpthis; + + char* buffer=*pbuffer; + + int i, j, size_tria=0; + int ig; + + int num_ints=0, num_doubles=0; + + + /* first marshall the embedded Tria object: ... */ + + TriaMarshall( &buffer, &this->tria, &size_tria); + + + /* ... then marshall the TriaElement member data: */ + + for ( i=0; i<3; i++) { + cellmar( &this->internalgrid_indices[i], &buffer, INTEGER_FLAG); + cellmar( &this->internalgrid_indicesb[i], &buffer, INTEGER_FLAG); + num_ints++; + } + /* marshall InternalGrid pointers even though pointer information won't be + valid after demarshalling (assumes pointer can fit into integer space): */ + for ( i=0; i<3; i++) { + cellmar( &this->pg[i], &buffer, INTEGER_FLAG); + num_ints++; + } + + cellmar( &this->property_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->property_enum, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->property, &buffer, INTEGER_FLAG); num_ints++; + + cellmar( &this->matice_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matice_enum, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matice, &buffer, INTEGER_FLAG); num_ints++; + + cellmar( &this->matpar_index, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matpar_enum, &buffer, INTEGER_FLAG); num_ints++; + cellmar( &this->matpar, &buffer, INTEGER_FLAG); num_ints++; + + *size = size_tria + num_ints*sizeof(int) + num_doubles*sizeof(double); + *pbuffer = buffer; + + return; +} + + +/*-------------------------------------------------- + TriaElementDemarshall + --------------------------------------------------*/ + +void* TriaElementDemarshall( DataSet* ds, char* *pbuffer, int machine_flag) +{ + int i, j, ig; + Tria* ptria=NULL; /* for embedded object */ + + char* buffer=*pbuffer; + TriaElement* this=NULL; + + + /* first demarshall the embedded Tria object: ... */ + + ptria = TriaDemarshall( NULL, &buffer, machine_flag); + /* create the derived object: */ + this = NewTriaElement( ptria); + /* and free tria now that it's embedded in the derived object: */ + DeleteTria( (void**)&ptria); + + + /* ... then demarshall the TriaElement member data: */ + + for ( i=0; i<3; i++) { + celldem( &this->internalgrid_indices[i], &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->internalgrid_indicesb[i], &buffer, INTEGER_FLAG, machine_flag); + } + + /* though they're meaningless, demarshall InternalGrid pointers: */ + for ( i=0; i<3; i++) { /* order+1 due to collapsed node formulation */ + celldem( &this->pg[i], &buffer, INTEGER_FLAG, machine_flag); + } + + celldem( &this->property_index, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->property_enum, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->property, &buffer, INTEGER_FLAG, machine_flag); + + + celldem( &this->property_index, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->property_enum, &buffer, INTEGER_FLAG, machine_flag); + celldem( &this->property, &buffer, INTEGER_FLAG, machine_flag); + + /*Set pointers to NULL :*/ + for (i=0;i<3;i++){ + this->pg[i]=NULL; + } + this->property=NULL; + this->matice=NULL; + this->matpar=NULL; + + if (ds) TriaElementAddToDataSet( ds, this); + *pbuffer=buffer; + + return (void*)this; +} + +/*-------------------------------------------------- + TriaElementGetID + --------------------------------------------------*/ + +int TriaElementGetID( void* vpthis ) +{ + /* vpthis for polymorphic function compatibility */ + TriaElement* this = (TriaElement*)vpthis; + + return TriaGetID( &this->tria); +} + + + +/*-------------------------------------------------- + TriaElementSetVirtualFunctions + --------------------------------------------------*/ + +int TriaElementSetVirtualFunctions( void* vp) +{ + + /* VirtualFunctions implemented as a "friend" class: */ + + /* void* for polymorphic function compatibility */ + VirtualFunctions* vf = (VirtualFunctions*)vp; + + if ( VirtualFunctionsInit(vf) ) { + + /* set TriaElement-specific virtual function pointers: */ + + /* general: */ + + vf->Configure = &TriaElementConfigure; + + /* grids, connectivity, basis conversions for element types: */ + + /* element matrix generation: */ + vf->CreateKMatrix = &TriaElementCreateKMatrix; + + /* element load generation: */ + vf->CreatePVector = &TriaElementCreatePVector; + + /* velocity difference generation: */ + vf->CreateDuVector = &TriaElementCreateDuVector; + vf->CreateGradjVectors = &TriaElementCreateGradjVectors; + + /* velocity misfit: */ + vf->Misfit=&TriaElementMisfit; + + + + return 1; + } + else + return 0; + +} + +/*-------------------------------------------------- + TriaElementConfigure + --------------------------------------------------*/ + +int TriaElementConfigure( void* vpthis, int num_datasets, ...) +{ + + /* set up data connections to all necessary, related objects: */ + + TriaElement* this = (TriaElement*)vpthis; + + int i,j, found, foundgrid, noerr=1, first=1; + + /* expected DataSet input: */ + + DataSet* pgpdt=NULL; + DataSet* pgpdtb=NULL; + DataSet* pmpt=NULL; + DataSet* pest=NULL; + + /* variable argument list handling: */ + + va_list args; + DataSet* *arglist=NULL; + + #ifdef _DEBUG_ + printf("virtual function: TriaElementConfigure\n"); + #endif + + + if (this) { + + if (!num_datasets) { + + printf(" %s error:\n", __func__); + printf(" null input for number of datasets;\n"); + printf(" could not establish grid, property and material data 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<3;i++){ + this->pg[i]=NULL; + } + + /* grid connections based on DataSet of type "gpdt": */ + + #ifndef _PARALLEL_ + found=0; + for ( i=0; itria); + + for ( i=0; i<3; 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) { + 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, TriaGetID(&this->tria), TriaGetName(&this->tria)); + } + } + } + } + + } + else { + + if (first) { + printf(" %s error:\n", __func__); + first=0; + } + printf(" could not find DataSet of type \"gpdt\", grid object\n"); + printf(" pointers not established for tria element id = %d\n", + TriaElementGetID(this)); + noerr = 0; + } + #else //ifdef _PARALLEL_ + + /*Configure the element when we are running on a parallel cluster: in this case, + we have partitioned the elements and the grids correspondingly. Every node has + its own est, 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 + element, 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; itria); + + for ( i=0; i<3; 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_indices[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, TriaGetID(&this->tria), TriaGetName(&this->tria)); + } + } + } + } + } + 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 tria element id = %d\n", + TriaElementGetID(this)); + noerr = 0; + } + + #endif + + /* material property connections based on DataSet of type "mpt": */ + + found=0; + for ( i=0; itria); + + /* as for grids, check to see if a nonzero matice index still + points to the matice entry of interest in mpt, otherwise do + a seek for it: */ + + if ( this->matice_index != UNDEF) { + + + if ( this->matice_enum == DataSetGetEnum( pmpt, this->matice_index) + && + current_matice_id == DataSetGetObjectID( pmpt, this->matice_index) ) { + + /* index is still valid; just get matice entry pointer directly: */ + + this->matice = DataSetGetObjectPtr( pmpt, this->matice_index); + } + } + + /* if we haven't yet "recovered" the matice object pointer, go seek it out: */ + + if ( NULL == this->matice) { + + noerr *= DataSetSeek( &this->matice, &this->matice_index, + pmpt, current_matice_id, 1, MaticeEnum() ); /* 0 to generalise later to Pcomp, etc. ... */ + + if (!noerr) { + if (first) { + printf(" %s error:\n", __func__); + first = 0; + } + printf(" could not find matice entry with id of %d for element %d of type \"%s\"\n", + current_matice_id, TriaGetID(&this->tria), TriaGetName(&this->tria)); + } + else { + + /* though the answer for now is just MaticeEnum(), for generality go ahead and + ask the question, "what type of matice has been referenced?" */ + + this->matice_enum = DataSetGetEnum( *(arglist+i), this->matice_index); + } + } + + int current_matpar_id = TriaGetMaterialParID(&this->tria); + + /* as for grids, check to see if a nonzero matpar index still + points to the matpar entry of interest in mpt, otherwise do + a seek for it: */ + + if ( this->matpar_index != UNDEF) { + + + if ( this->matpar_enum == DataSetGetEnum( pmpt, this->matpar_index) + && + current_matpar_id == DataSetGetObjectID( pmpt, this->matpar_index) ) { + + /* index is still valid; just get matpar entry pointer directly: */ + + this->matpar = DataSetGetObjectPtr( pmpt, this->matpar_index); + } + } + + /* if we haven't yet "recovered" the matpar object pointer, go seek it out: */ + + if ( NULL == this->matpar) { + + noerr *= DataSetSeek( &this->matpar, &this->matpar_index, + pmpt, current_matpar_id, 1, MatparEnum() ); /* 0 to generalise later to Pcomp, etc. ... */ + + if (!noerr) { + if (first) { + printf(" %s error:\n", __func__); + first = 0; + } + printf(" could not find matpar entry with id of %d for element %d of type \"%s\"\n", + current_matpar_id, TriaGetID(&this->tria), TriaGetName(&this->tria)); + } + else { + + /* though the answer for now is just MatparEnum(), for generality go ahead and + ask the question, "what type of matpar has been referenced?" */ + + this->matpar_enum = DataSetGetEnum( *(arglist+i), this->matpar_index); + } + } + } + else { + + if (first) { + printf(" %s error:\n", __func__); + first=0; + } + printf(" could not find DataSet of type \"mpt\"; matice and matpar relationships\n"); + printf(" not established for tria element id = %d\n", + TriaElementGetID(this)); + noerr = 0; + } + + + #ifdef _DEBUG_ + printf(" element data after all lookups:\n"); + TriaElementEcho( this); + #endif + + + xfree((void**)&arglist); + } + else { + + /* quiet return */ + ; + } + cleanup_and_return: + + return noerr; + +} + +/*-------------------------------------------------- + TriaElementAddToDataSet + --------------------------------------------------*/ + +int TriaElementAddToDataSet( DataSet* dataset, TriaElement* 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, TriaElementEnum(), (void*)this, + &DeleteTriaElement, &TriaElementEcho, &TriaElementSizeOfInDoubles, + &TriaElementShrinkWrap, &TriaElementUnwrap, + &TriaElementMarshall, + &TriaElementGetID, &TriaElementSetVirtualFunctions); + + return 1; +} + +/*-------------------------------------------------- + TriaElementCreateKMatrix + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreateKMatrix" +int TriaElementCreateKMatrix( ElemMatrix* *pKe_gg, ElemMatrix** pKTe_gg, void* vpthis, ParameterInputs* inputs, int K_flag, int KT_flag, int analysis_type){ + + + int noerr=1; + /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ + if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ + + noerr=TriaElementCreateKMatrixDiagnosticHoriz( pKe_gg, vpthis, inputs); + + } + else{ + printf("%s%s%i%s\n",__FUNCT__," error message: analysis ",analysis_type," not supported yet"); + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + TriaElementCreatePVector + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreatePVector" + +int TriaElementCreatePVector( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs, int analysis_type){ + + int noerr=1; + + /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ + if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){ + noerr=TriaElementCreatePVectorDiagnosticHoriz( ppe_g,vpthis,inputs); + } + else{ + printf("%s%s%i%s\n",__FUNCT__," error message: analysis ",analysis_type," not supported yet"); + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + TriaElementCreateKMatrixDiagnosticHoriz + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreateKMatrixDiagnosticHoriz" +int TriaElementCreateKMatrixDiagnosticHoriz( ElemMatrix* *pKe_gg, void* vpthis, ParameterInputs* inputs){ + + + int i,j,il,counter; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + TriaElement* this = NULL; + Matice * matice = NULL; + Matpar * matpar = NULL; + + /* local declarations */ + + /* output: */ + ElemMatrix* Ke_gg = NULL; + + + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + + /* specific gaussian point: */ + double gauss_weight; + double gauss_l1l2l3[3]; + + /*drag: */ + double pcoeff,qcoeff; + double rcoeff,scoeff; + double velocity_x,velocity_y,velocity_mag; + + /* material data: */ + double viscosity; //viscosity + double old_viscosity; + double viscosity_change; + double B_param; + double gravity=9.8; + double rho_ice,rho_water; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + double old_epsilon[3]; + + /* matrices: */ + double B[3][NDOF2*numgrids]; + double Bprime[3][NDOF2*numgrids]; + double L[2][NDOF2*numgrids]; + double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. + double D_scalar; + + double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + double DL_scalar; + + double Ke_gg_gaussian[NDOF2*numgrids][NDOF2*numgrids]; //stiffness matrix evaluated at the gaussian point. + double Ke_gg_drag_gaussian[NDOF2*numgrids][NDOF2*numgrids]; //stiffness matrix contribution from drag + double Jdet; + + /*slope: */ + double slope[NDOF2]={0.0,0.0}; + double slope_magnitude; + + /*input parameters for structural analysis (diagnostic): */ + double* thickness_param=NULL; + double thickness_list[numgrids]; + double thickness; + double* surface_param=NULL; + double surface_list[numgrids]; + double* bed_param=NULL; + double bed_list[numgrids]; + double* flow_law_param=NULL; + double B_list[numgrids]; + double* basal_drag=NULL; + double K_list[numgrids]; + double* velocity=NULL; + double vxvy_list[numgrids][2]; + double* temperature_average_param=NULL; + double temperature_average_list[numgrids]; + double temperature_average; + + /*friction: */ + double alpha2_list[numgrids]={0.0,0.0,0.0}; + double alpha2; + + double MAXSLOPE=.06; // 6 % + double MOUNTAINKEXPONENT=10; + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids + * per element: */ + numdof=numgrids*NDOF2; + + /*Some pointer intialization: */ + this = (TriaElement*)vpthis; + matice=this->matice; + matpar=this->matpar; + + /*First check that acceleration is not on: */ + if (this->tria.acceleration){ + *pKe_gg=NULL; + return; + } + + /* Allocate new stiffness element matrix: */ + Ke_gg=NewElemMatrix(numdof); + + /*Recover material: */ + gravity=matpar->g; + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + + /*recover extra inputs from users, at current convergence iteration: */ + basal_drag=ParameterInputsRecover(inputs,"drag"); + thickness_param=ParameterInputsRecover(inputs,"thickness"); + surface_param=ParameterInputsRecover(inputs,"surface"); + bed_param=ParameterInputsRecover(inputs,"bed"); + temperature_average_param=ParameterInputsRecover(inputs,"temperature_average"); + flow_law_param=ParameterInputsRecover(inputs,"B"); + velocity=ParameterInputsRecover(inputs,"velocity"); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not get element grid data!"); + goto cleanup_and_return; + } + + /* Build the row index vector (same as column index vector). + * At the same time, get element exterior data for each grid: */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=dof; + if(velocity){ + vxvy_list[i][j]=*(velocity+dof); + } + else{ + vxvy_list[i][j]=0; + } + } + dof=structural_dof_list[0]; + if(flow_law_param) + B_list[i] = *(flow_law_param+dof); + else + B_list[i] =0; + + if(basal_drag){ + K_list[i]=*(basal_drag+dof); + } + else{ + K_list[i]=this->tria.k[i]; + } + if(thickness_param){ + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->tria.h[i]; + } + if(surface_param){ + surface_list[i]= *(surface_param+dof); + } + else{ + surface_list[i]= this->tria.s[i]; + } + if(bed_param){ + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]= this->tria.b[i]; + } + + if(temperature_average_param) temperature_average_list[i]=*(temperature_average_param+dof); + } + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf("El id %i Rank %i TriaElemnet input list before gaussian loop: \n",ELID,RANK); + printf(" rho_ice: %g \n",rho_ice); + printf(" gravity: %g \n",gravity); + printf(" rho_water: %g \n",rho_water); + printf(" Velocity: \n"); + for (i=0;itria.shelf && (this->tria.friction_type==2)){ + + /*Allocate friction object: */ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=gravity; + friction->rho_ice=rho_ice; + friction->rho_water=rho_water; + friction->K=&K_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=this->tria.p; + friction->q=this->tria.q; + + /*Compute alpha2_list: */ + FrictionGetAlpha2(&alpha2_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + } + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf(" alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); + } + #endif + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf(" gaussian points: \n"); + for(i=0;i 6% for this element, it means we are on a mountain. In this particular case, + //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ + if(!this->tria.shelf){ + TriaElementGetParameterDerivativeValue(&slope[0], &surface_list[0],&xyz_list[0][0], gauss_l1l2l3); + slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); + + if (slope_magnitude>MAXSLOPE){ + alpha2_list[0]=pow(10,MOUNTAINKEXPONENT); + alpha2_list[1]=pow(10,MOUNTAINKEXPONENT); + alpha2_list[2]=pow(10,MOUNTAINKEXPONENT); + } + } + + /*Get strain rate from velocity: */ + TriaElementGetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); + + //Update material if temperature is provided. + if(temperature_average_param){ + TriaElementGetParameterValue(&temperature_average, &temperature_average_list[0],gauss_l1l2l3); + B_param=Paterson(temperature_average); + MaticeSetFlowLawParameter(this->matice,B_param); + } + //Update material if B is provided. + if(flow_law_param){ + TriaElementGetParameterValue(&B_param, &B_list[0],gauss_l1l2l3); + MaticeSetFlowLawParameter(this->matice,B_param); + } + + /*Get viscosity: */ + MaticeGetViscosity2d(&viscosity, this->matice, &epsilon[0]); + + + /* Get Jacobian determinant: */ + TriaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + onto this scalar matrix, so that we win some computational time: */ + D_scalar=viscosity*thickness*gauss_weight*Jdet; + for (i=0;i<3;i++){ + D[i][i]=D_scalar; + } + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf(" gaussian loop %i\n",ig); + printf(" thickness %g\n",thickness); + printf(" slope [%g,%g]\n",slope[0],slope[1]); + printf(" alpha2_list [%g,%g,%g]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]); + printf(" epsilon [%g,%g,%g]\n",epsilon[0],epsilon[1],epsilon[2]); + printf(" Matice: \n"); + MaticeEcho(this->matice); + printf("\n viscosity: %g \n",viscosity); + printf(" jacobian: %g \n",Jdet); + printf(" gauss_weight: %g \n",gauss_weight); + } + #endif + + /*Get B and Bprime matrices: */ + TriaElementGetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3); + TriaElementGetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3); + + /*Get L matrix if viscous basal drag present: */ + if((this->tria.friction_type==2) && (!this->tria.shelf)){ + TriaElementGetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF2); + } + /* Do the triple product tB*D*Bprime: */ + TripleMultiply( &B[0][0],3,numdof,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ + for( i=0; inrows; i++){ + for (j=0;jnrows;j++){ + *(Ke_gg->terms+Ke_gg->nrows*i+j)+=Ke_gg_gaussian[i][j]; + } + } + + /*Now, take care of the basal friction if there is any: */ + if((!this->tria.shelf) && (this->tria.friction_type==2)){ + + TriaElementGetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3); + + if (velocity){ + DL_scalar=alpha2*gauss_weight*Jdet; + } + else{ + DL_scalar=0; + } + + for (i=0;i<2;i++){ + DL[i][i]=DL_scalar; + } + + /* Do the triple producte tL*D*L: */ + TripleMultiply( &L[0][0],2,numdof,1, + &DL[0][0],2,2,0, + &L[0][0],2,numdof,0, + &Ke_gg_drag_gaussian[0][0],0); + + for( i=0; inrows; i++){ + for (j=0;jnrows;j++){ + *(Ke_gg->terms+Ke_gg->nrows*i+j)+=Ke_gg_drag_gaussian[i][j]; + } + } + } + + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf(" alpha2 %g\n",alpha2); + printf(" B:\n"); + for(i=0;i<3;i++){ + for(j=0;jterms:\n"); + for( i=0; inrows; i++){ + for (j=0;jnrows;j++){ + printf("%g ",*(Ke_gg->terms+Ke_gg->nrows*i+j)); + } + printf("\n"); + } + printf(" Ke_gg->row_indices:\n"); + for( i=0; inrows; i++){ + printf("%i ",Ke_gg->row_indices[i]); + } + + } + #endif + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build stiffness matrix."); + return 0; + } + else{ + /*Assign output pointer: */ + *pKe_gg=Ke_gg; + + return 1; + } + +} + +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreatePVectorDiagnosticHoriz" + +int TriaElementCreatePVectorDiagnosticHoriz( ElemVector* *ppe_g, void* vpthis, ParameterInputs* inputs){ + + int i,j; + int noerr=1; + + /* vpthis for polymorphic function compatibility */ + TriaElement* this = NULL; + Matice * matice = NULL; + Matpar * matpar = NULL; + + /* output: */ + ElemVector* pe_g = NULL; + + int numdof; + int* structural_dof_list = NULL; + int dof; + + + /* parameters: */ + double rho_ice; + double gravity=9.8; + double plastic_stress; + double slope[NDOF2]; + double driving_stress_baseline; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /*element vector at the gaussian points: */ + double pe_g_gaussian[numgrids*NDOF2]; + + /*input parameters for structural analysis (diagnostic): */ + double* thickness_param=NULL; + double thickness_list[numgrids]; + double thickness; + double* surface_param=NULL; + double surface_list[numgrids]; + double* bed_param=NULL; + double bed_list[numgrids]; + double* basal_drag_param=NULL; + double K_list[numgrids]; + + + /*Some pointer intialization: */ + this = (TriaElement*)vpthis; + matice=this->matice; + matpar=this->matpar; + + /*First check that acceleration is not on: */ + if (this->tria.acceleration){ + *ppe_g=NULL; + return; + } + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids + * per element: */ + numdof=numgrids*NDOF2; + + /* Allocate new load element vector: */ + pe_g=NewElemVector(numdof); + + /* recover material parameters: */ + gravity=matpar->g; + rho_ice=matpar->rho_ice; + + /* recover input parameters: */ + thickness_param=ParameterInputsRecover(inputs,"thickness"); + surface_param=ParameterInputsRecover(inputs,"surface"); + bed_param=ParameterInputsRecover(inputs,"bed"); + basal_drag_param=ParameterInputsRecover(inputs,"drag"); + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], &cid_list[0], this->pg, numgrids, 1 ); + + /* Build the row index vector : */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=structural_dof_list[j]; + } + if(thickness_param){ + dof=structural_dof_list[0]; + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->tria.h[i]; + } + + if(surface_param){ + dof=structural_dof_list[0]; + surface_list[i]=*(surface_param+dof); + } + else{ + surface_list[i]=this->tria.s[i]; + } + if(bed_param){ + dof=structural_dof_list[0]; + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]=this->tria.b[i]; + } + if(basal_drag_param){ + dof=structural_dof_list[0]; + K_list[i]=*(basal_drag_param+dof); + } + else{ + K_list[i]=this->tria.k[i]; + } + } + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf("gravity %g\n",gravity); + printf("rho_ice %g\n",rho_ice); + printf("thickness_list [%g,%g,%g]\n",thickness_list[0],thickness_list[1],thickness_list[2]); + printf("surface_list [%g,%g,%g]\n",surface_list[0],surface_list[1],surface_list[2]); + printf("bed_list [%g,%g,%g]\n",bed_list[0],bed_list[1],bed_list[2]); + printf("K_list [%g,%g,%g]\n",K_list[0],K_list[1],K_list[2]); + } + #endif + + + /* Get gaussian points and weights: */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf(" gaussian points: \n"); + for(i=0;itria.friction_type==1){ + TriaElementGetParameterValue(&plastic_stress, &K_list[0],gauss_l1l2l3); + } + + /* Get Jacobian determinant: */ + TriaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + + /*Get nodal functions: */ + TriaElementGetNodalFunctions(l1l2l3, gauss_l1l2l3); + + /*Compute driving stress: */ + driving_stress_baseline=rho_ice*gravity*thickness; + + + #ifdef _DEBUGELEMENTS_ + if(my_rank==RANK && TriaElementGetID(this)==ELID){ + printf(" gaussian %i\n",ig); + printf(" thickness %g\n",thickness); + printf(" slope(%g,%g)\n",slope[0],slope[1]); + printf(" Jdet %g\n",Jdet); + printf(" gaussweigth %g\n",gauss_weight); + printf(" l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]); + if(this->tria.friction_type==1)printf(" plastic_stress(%g)\n",plastic_stress); + } + #endif + + /*Build pe_g_gaussian vector: */ + if(this->tria.friction_type==1){ + for (i=0;inrows; i++){ + *(pe_g->terms+i)+=*(pe_g_gaussian+i); + } + } //for (ig=0; igterms\n",ig); + for( i=0; inrows; i++){ + printf("%g ",*(pe_g->terms+i)); + } + printf("\n"); + printf(" pe_g->row_indices\n",ig); + for( i=0; inrows; i++){ + printf("%i ",*(pe_g->row_indices+i)); + } + } + #endif + + + cleanup_and_return: + + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build load vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(pe_g); + #endif + + /*Assign output pointer: */ + *ppe_g=pe_g; + return 1; + } +} + + + +/*-------------------------------------------------- + TriaElementGetParameterValue + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementGetParameterValue" +int TriaElementGetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){ + + /*From grid values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian + * point specifie by gauss_l1l2l3: */ + + int noerr=1; + + /*nodal functions: */ + double l1l2l3[3]; + + /*output: */ + double p; + + noerr = TriaElementGetNodalFunctions(l1l2l3, gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2]; + + + cleanup_and_return: + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not return parameter value."); + return 0; + } + else{ + *pp=p; + return 1; + } +} + + +/*-------------------------------------------------- + TriaElementGetParameterDerivativeValue + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementGetParameterDerivativeValue" +int TriaElementGetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3){ + + /*From grid values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian + * point specified by gauss_l1l2l3: + * dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx + * dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx + * + * p is a vector of size 2x1 already allocated. + */ + + int noerr=1; + + double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system. + + /*Get dh1dh2dh3 in basic coordinate system: */ + noerr=TriaElementGetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + *(p+0)=plist[0]*dh1dh2dh3_basic[0][0]+plist[1]*dh1dh2dh3_basic[0][1]+plist[2]*dh1dh2dh3_basic[0][2]; + *(p+1)=plist[0]*dh1dh2dh3_basic[1][0]+plist[1]*dh1dh2dh3_basic[1][1]+plist[2]*dh1dh2dh3_basic[1][2]; + + + cleanup_and_return: + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not return parameter derivative vector."); + return 0; + } + else{ + return 1; + } +} + +/*-------------------------------------------------- + TriaElementGetStrainRate + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementGetStrainRate" +int TriaElementGetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3){ + + int noerr=1; + int i; + + double B[3][NDOF2*numgrids]; + + /*Get B matrix: */ + noerr=TriaElementGetB(&B[0][0], xyz_list, gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + #ifdef _DEBUG_ + printf("B for grid1 : [ %lf %lf \n",B[0][0],B[0][1]); + printf(" [ %lf %lf \n",B[1][0],B[1][1]); + printf(" [ %lf %lf ]\n",B[2][0],B[2][1]); + + printf("B for grid2 : [ %lf %lf \n",B[0][2],B[0][3]); + printf(" [ %lf %lf \n",B[1][2],B[1][3]); + printf(" [ %lf %lf ]\n",B[2][2],B[2][3]); + + printf("B for grid3 : [ %lf %lf \n",B[0][4],B[0][5]); + printf(" [ %lf %lf \n",B[1][4],B[1][5]); + printf(" [ %lf %lf ]\n",B[2][4],B[2][5]); + + for (i=0;itria.fit; + } + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + + /* Build the row index vector : */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i*NDOF2+j)=dof; + } + dof=structural_dof_list[0]; + vx_list[i]=*(velocity+dof); + obs_vx_list[i]=*(obs_velocity+dof); + dof=structural_dof_list[1]; + vy_list[i]=*(velocity+dof); + obs_vy_list[i]=*(obs_velocity+dof); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + #ifdef _DEBUG_ + for (i=0;itria.meanvel; + epsvel=this->tria.epsvel; + scalex=pow(meanvel/(obs_velocity_x+epsvel),2); + scaley=pow(meanvel/(obs_velocity_y+epsvel),2); + if(obs_velocity_x==0)scalex=0; + if(obs_velocity_y==0)scaley=0; + for (i=0;itria.meanvel; + epsvel=this->tria.epsvel; + velocity_mag=sqrt(pow(velocity_x,2)+pow(velocity_y,2))+epsvel; //epsvel to avoid velocity being nil. + obs_velocity_mag=sqrt(pow(obs_velocity_x,2)+pow(obs_velocity_y,2))+epsvel; //epsvel to avoid observed velocity being nil. + scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + for (i=0;inrows; i++){ + *(due_g->terms+i)+=*(due_g_gaussian+i); + } + } + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build velocity difference vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(due_g); + #endif + + /*Assign output pointer: */ + *pdue_g=due_g; + + return 1; + } +} + + +/*-------------------------------------------------- + TriaElementCreateGradjVectors + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreateGradVectors" +int TriaElementCreateGradjVectors( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type,char* control_type){ + + int noerr=1; + + if strcasecmp_eq(control_type,"drag"){ + noerr=TriaElementCreateGradjVectorsDrag( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + } + else if strcasecmp_eq(control_type,"B"){ + noerr=TriaElementCreateGradjVectorsB( pgradje_g, vpthis, velocity, adjoint, inputs,analysis_type); + } + else{ + printf("%s%s\n",__FUNCT__," error message: control type not yet supported!"); + noerr=0; + } + return noerr; +} + +/*-------------------------------------------------- + TriaElementCreateGradjVectorsDrag + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreateGradVectorsDrag" +int TriaElementCreateGradjVectorsDrag( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ + + int i,j; + int noerr; + + /* vpthis for polymorphic function compatibility */ + TriaElement* this = NULL; + Matice * matice = NULL; + Matpar * matpar = NULL; + + /* local declarations */ + + /* output: */ + ElemVector* gradje_g = NULL; + + int num_dof; + int* structural_dof_list = NULL; + int dof; + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + double vx_list[numgrids]; + double vy_list[numgrids]; + double vxvy_list[numgrids][2]; + double adjx_list[numgrids]; + double adjy_list[numgrids]; + + double* basal_drag_param=NULL; + double K_list[numgrids]; + double k; + + double* bed_param=NULL; + double bed_list[numgrids]; + double* thickness_param=NULL; + double thickness_list[numgrids]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + + /* specific gaussian point: */ + double gauss_weight; + double gauss_l1l2l3[3]; + + /* parameters: */ + double viscosity2; + double dvx[NDOF2]; + double dvy[NDOF2]; + double dadjx[NDOF2]; + double dadjy[NDOF2]; + double vx,vy; + double lambda,mu; + double bed,thickness,Neff; + + /*drag: */ + double pcoeff,qcoeff; + double gravity=9.8; + double rho_ice,rho_water; + double alpha_complement_list[numgrids]; + double alpha_complement; + + + /*element vector at the gaussian points: */ + double gradje_g_gaussian[numgrids]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + + /*Some pointer intialization: */ + this = (TriaElement*)vpthis; + matice=this->matice; + matpar=this->matpar; + + /*recover material parameters: */ + rho_ice=matpar->rho_ice; + rho_water=matpar->rho_water; + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids + * per element: */ + num_dof=numgrids*NDOF2; + + /* Allocate element vector: */ + gradje_g=NewElemVector(numgrids); + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + + /* recover input parameters: */ + basal_drag_param=ParameterInputsRecover(inputs,"drag"); + bed_param=ParameterInputsRecover(inputs,"bed"); + thickness_param=ParameterInputsRecover(inputs,"thickness"); + + /* Build the row index vector : */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i)=dof; + + if(basal_drag_param){ + K_list[i]=*(basal_drag_param+dof); + } + else{ + K_list[i]=this->tria.k[i]; + } + if(bed_param){ + bed_list[i]=*(bed_param+dof); + } + else{ + bed_list[i]=this->tria.b[i]; + } + if(thickness_param){ + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->tria.h[i]; + } + + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); + #ifdef _DEBUG_ + for (i=0;itria.shelf && (this->tria.friction_type==2)){ + + /*Allocate friction object: */ + Friction* friction=NewFriction(); + + /*Initialize all fields: */ + friction->element_type=xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=gravity; + friction->rho_ice=rho_ice; + friction->rho_water=rho_water; + friction->K=&K_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->velocities=&vxvy_list[0][0]; + friction->p=this->tria.p; + friction->q=this->tria.q; + + if(friction->p!=1){ + printf("Non-linear friction not supported yet in control methods!\n"); + noerr=0;goto cleanup_and_return; + } + + /*Compute alpha complement: */ + FrictionGetAlphaComplement(&alpha_complement_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + } + else{ + alpha_complement_list[0]=0; + alpha_complement_list[1]=0; + alpha_complement_list[2]=0; + } + + /*Recover alpha_complement and k: */ + noerr = TriaElementGetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3); + noerr = TriaElementGetParameterValue(&k, &K_list[0],gauss_l1l2l3); + #ifdef _DEBUG_ + printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,k); + #endif + + /*recover lambda and mu: */ + noerr = TriaElementGetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + noerr = TriaElementGetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + #ifdef _DEBUG_ + printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu); + #endif + + /*recover vx and vy: */ + noerr = TriaElementGetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + noerr = TriaElementGetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + #ifdef _DEBUG_ + printf("Velocity vector %20.20lf %20.20lf\n",vx,vy); + #endif + + /* Get Jacobian determinant: */ + noerr=TriaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + #ifdef _DEBUG_ + printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet); + #endif + + /* Get nodal functions value at gaussian point:*/ + noerr = TriaElementGetNodalFunctions(l1l2l3, gauss_l1l2l3); + if(!noerr)goto cleanup_and_return; + + + /*Build gradje_g_gaussian vector: */ + for (i=0;inrows; i++){ + *(gradje_g->terms+i)+=*(gradje_g_gaussian+i); + } + } + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build velocity difference vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(gradje_g); + #endif + + + /*Assign output pointer: */ + *pgradje_g=gradje_g; + + return 1; + } +} + + +/*-------------------------------------------------- + TriaElementCreateGradjVectorsB + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementCreateGradjVectorsB" +int TriaElementCreateGradjVectorsB( ElemVector* *pgradje_g, void* vpthis, double* velocity, double* adjoint, ParameterInputs* inputs,int analysis_type){ + + int i,j; + int noerr; + + /* vpthis for polymorphic function compatibility */ + TriaElement* this = NULL; + + /* local declarations */ + + /* output: */ + ElemVector* gradje_g = NULL; + + int num_dof; + int* structural_dof_list = NULL; + int dof; + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + double vx_list[numgrids]; + double vy_list[numgrids]; + double vxvy_list[numgrids][2]; + double adjx_list[numgrids]; + double adjy_list[numgrids]; + + double* thickness_param=NULL; + double thickness_list[numgrids]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + + /* specific gaussian point: */ + double gauss_weight; + double gauss_l1l2l3[3]; + + /* parameters: */ + double viscosity2; + double dvx[NDOF2]; + double dvy[NDOF2]; + double dadjx[NDOF2]; + double dadjy[NDOF2]; + double vx,vy; + double lambda,mu; + double thickness; + + + /*element vector at the gaussian points: */ + double gradje_g_gaussian[numgrids]; + + /* Jacobian: */ + double Jdet; + + /*nodal functions: */ + double l1l2l3[3]; + + /* strain rate: */ + double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + /*Some pointer intialization: */ + this = (TriaElement*)vpthis; + + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids + * per element: */ + num_dof=numgrids*NDOF2; + + /* Allocate element vector: */ + gradje_g=NewElemVector(numgrids); + + /* Get all element grid data */ + noerr = GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); if(!noerr)goto cleanup_and_return; + + /* recover input parameters: */ + thickness_param=ParameterInputsRecover(inputs,"drag"); + + /* Build the row index vector : */ + for (i=0;ipg[i]); + for (j=0;jrow_indices+i)=dof; + + if(thickness_param){ + thickness_list[i]=*(thickness_param+dof); + } + else{ + thickness_list[i]=this->tria.h[i]; + } + + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); + #ifdef _DEBUG_ + for (i=0;ictrice3.shelf){ + thickness=0; //this will ensure gradjbe_g_gaussian is nil; + }*/ + + /*Get strain rate, if velocity has been supplied: */ + noerr=TriaElementGetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + /*Get viscosity2: */ + noerr=MaticeGetViscosity2(&viscosity2, this->matice, &epsilon[0]); if(!noerr)goto cleanup_and_return; + + /*Get dvx: */ + noerr=TriaElementGetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + /*Get dvy: */ + noerr=TriaElementGetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + /*Get dadjx: */ + noerr=TriaElementGetParameterDerivativeValue(&dadjx[0], &adjx_list[0],&xyz_list[0][0], gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + /*Get dadjy: */ + noerr=TriaElementGetParameterDerivativeValue(&dadjy[0], &adjy_list[0],&xyz_list[0][0], gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + + /* Get Jacobian determinant: */ + noerr=TriaElementGetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + /* Get nodal functions value at gaussian point:*/ + noerr = TriaElementGetNodalFunctions(l1l2l3, gauss_l1l2l3); if(!noerr)goto cleanup_and_return; + + #ifdef _DEBUG_ + _printf_("viscosity2 %g thickness %g dvx [%g %g] dvy [%g %g] dadjx [%g %g] dadjy[%g %g]\n",viscosity2,thickness,dvx[0],dvx[1],dvy[0],dvy[1],dadjx[0],dadjx[1],dadjy[0],dadjy[1]); + #endif + + /*Build gradje_g_gaussian vector: */ + for (i=0;inrows; i++){ + *(gradje_g->terms+i)+=*(gradje_g_gaussian+i); + } + } + + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not build gradjB vector."); + return 0; + } + else{ + #ifdef _DEBUG_ + ElemVectorEcho(gradje_g); + #endif + /*Assign output pointer: */ + *pgradje_g=gradje_g; + + return 1; + } +} + + + + +/*-------------------------------------------------- + TriaElementMisfit + --------------------------------------------------*/ +#undef __FUNCT__ +#define __FUNCT__ "TriaElementMisfit" +int TriaElementMisfit( double* pJelem,void* vpthis, double* velocity, double* obs_velocity, ParameterInputs* inputs, int analysis_type){ + + int i,j; + int noerr; + + /* vpthis for polymorphic function compatibility */ + TriaElement* this = NULL; + + /* local declarations */ + + /* output: */ + double Jelem=0; + + int num_dof; + int* structural_dof_list = NULL; + int dof; + + /* grid data: */ + int cid_list[numgrids]; + double xyz_list[numgrids][3]; + double T_bg_list[numgrids][9]; + double vx_list[numgrids]; + double vy_list[numgrids]; + double obs_vx_list[numgrids]; + double obs_vy_list[numgrids]; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + + /* specific gaussian point: */ + double gauss_weight; + double gauss_l1l2l3[3]; + + /* parameters: */ + double velocity_x,velocity_y,velocity_mag; + double obs_velocity_x,obs_velocity_y,obs_velocity_mag; + + /* Jacobian: */ + double Jdet; + + /*relative and logarithmic control method :*/ + double meanvel=0; + double epsvel=0; + double scalex=1; + double scaley=1; + double* fit_param=NULL; + double fit=-1; + + + /*Some pointer intialization: */ + this = (TriaElement*)vpthis; + + /* The number of g-set dof for this element is 2*3 (vx and vy, horizontal velocities for ice, and 3 grids + * per element: */ + num_dof=numgrids*NDOF2; + + /* Get all element grid data */ + noerr *= GetElementGridData( &xyz_list[0][0], &T_bg_list[0][0], cid_list, this->pg, numgrids, 1 ); + + /* Recover input data: */ + fit_param=ParameterInputsRecover(inputs,"fit"); + if(fit_param){ + fit=*fit_param; + } + else{ + fit=this->tria.fit; + } + + /* Build the row index vector : */ + for (i=0;ipg[i]); + dof=structural_dof_list[0]; + vx_list[i]=*(velocity+dof); + obs_vx_list[i]=*(obs_velocity+dof); + dof=structural_dof_list[1]; + vy_list[i]=*(velocity+dof); + obs_vy_list[i]=*(obs_velocity+dof); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + #ifdef _DEBUG_ + for (i=0;itria.meanvel; + epsvel=this->tria.epsvel; + scalex=pow(meanvel/(obs_velocity_x+epsvel),2); + scaley=pow(meanvel/(obs_velocity_y+epsvel),2); + if(obs_velocity_x==0)scalex=0; + if(obs_velocity_y==0)scaley=0; + Jelem+=.5*(scalex*pow((velocity_x-obs_velocity_x),2)+scaley*pow((velocity_y-obs_velocity_y),2))*Jdet*gauss_weight; + } + else if(fit==2){ + /*Logarithmic misfit: */ + meanvel=this->tria.meanvel; + epsvel=this->tria.epsvel; + velocity_mag=sqrt(pow(velocity_x,2)+pow(velocity_y,2))+epsvel; //epsvel to avoid velocity being nil. + obs_velocity_mag=sqrt(pow(obs_velocity_x,2)+pow(obs_velocity_y,2))+epsvel; //epsvel to avoid observed velocity being nil. + Jelem+=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2); + } + else{ + printf("%s%s\n",__FUNCT__," error message: unsupported type of fitting."); + noerr=0;goto cleanup_and_return; + } + + } + cleanup_and_return: + xfree((void**)&first_gauss_area_coord); + xfree((void**)&second_gauss_area_coord); + xfree((void**)&third_gauss_area_coord); + xfree((void**)&gauss_weights); + + if(!noerr){ + printf("%s%s\n",__FUNCT__," error message: could not compute velocity misfit."); + return 0; + } + else{ + /*Assign output pointer: */ + *pJelem=Jelem; + + return 1; + } +} +#undef NDOF1 +#undef NDOF2 +#undef numgrids + + +