Context Navigation

source: issm/trunk/src/c/objects/Pengrid.cpp@ 3621

Last change on this file since 3621 was 3621, checked in by Eric.Larour, 15 years ago
Almost there
File size: 22.6 KB

Rev	Line
[387]	1	/*!\file Pengrid.c
	2	* \brief: implementation of the Pengrid object
	3	*/
	4
	5
	6	#ifdef HAVE_CONFIG_H
	7	#include "config.h"
	8	#else
	9	#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
	10	#endif
	11
	12	#include "stdio.h"
	13	#include "./Pengrid.h"
	14	#include <string.h>
	15	#include "../EnumDefinitions/EnumDefinitions.h"
[3454]	16	#include "../include/macros.h"
[387]	17	#include "../shared/shared.h"
	18	#include "../include/typedefs.h"
[3332]	19	#include "../include/macros.h"
[3621]	20	#include "../DataSet/DataSet.h"
	21	#include "../DataSet/Inputs.h"
	22
[2907]	23	/Object constructors and destructor/
	24	/FUNCTION Pengrid::constructor {{{1/
[387]	25	Pengrid::Pengrid(){
[3621]	26	this->inputs=NULL;
	27	this->parameters=NULL;
	28
	29	/not active, not zigzagging: /
	30	active=0;
	31	zigzag_counter=0;
	32
[387]	33	}
[2907]	34	/}}}1/
[3621]	35	/FUNCTION Pengrid::Pengrid(int id, int node_ids int matpar_id){{{1/
	36	Pengrid::Pengrid(int pengrid_id,int pengrid_node_id, int pengrid_matpar_id):
	37	hnode(pengrid_node_ids,1),
	38	hmatice(&pengrid_matice_id,1),
	39	hmatpar(&pengrid_matpar_id,1)
	40	{
	41
	42	/all the initialization has been done by the initializer, just fill in the id: /
	43	this->id=pengrid_id;
	44	this->parameters=NULL;
	45	this->inputs=new Inputs();
	46
	47	/not active, not zigzagging: /
	48	active=0;
	49	zigzag_counter=0;
	50
	51	}
	52	/}}}/
	53	/FUNCTION Pengrid::Pengrid(int id, Hook hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* pengrid_inputs) {{{1*/
	54	Pengrid::Pengrid(int pengrid_id,Hook* pengrid_hnode, Hook* pengrid_hmatpar, Parameters* pengrid_parameters, Inputs* pengrid_inputs):
	55	hnode(pengrid_hnode),
	56	hmatpar(pengrid_hmatpar)
	57	{
	58
	59	/all the initialization has been done by the initializer, just fill in the id: /
	60	this->id=pengrid_id;
	61	if(pengrid_inputs){
	62	this->inputs=(Inputs*)pengrid_inputs->Copy();
	63	}
	64	else{
	65	this->inputs=new Inputs();
	66	}
	67	/point parameters: /
	68	this->parameters=pengrid_parameters;
	69
	70	/not active, not zigzagging: /
	71	active=0;
	72	zigzag_counter=0;
	73
	74	}
	75	/}}}/
[2212]	76	Pengrid::Pengrid(int pengrid_id, int pengrid_node_id,int pengrid_mparid, int pengrid_dof, int pengrid_active, double pengrid_penalty_offset,int pengrid_thermal_steadystate,int pengrid_stabilize_constraints){
[387]	77
	78	id=pengrid_id;
[495]	79	node_id=pengrid_node_id;
[483]	80	mparid=pengrid_mparid;
[387]	81	dof=pengrid_dof;
	82	active=pengrid_active;
	83	penalty_offset =pengrid_penalty_offset;
	84	thermal_steadystate=pengrid_thermal_steadystate;
[2212]	85	stabilize_constraints=pengrid_stabilize_constraints;
[387]	86
	87	node_offset=UNDEF;
	88	node=NULL;
[483]	89	matpar=NULL;
	90	matpar_offset=UNDEF;
[387]	91
[2212]	92
[387]	93	return;
	94	}
[2907]	95	/}}}1/
[3621]	96	/FUNCTION Pengrid::Pengrid(int i, IoModel iomodel){{{1*/
	97	Pengrid::Pengrid(int index, IoModel* iomodel){ //i is the element index
	98
	99	int i,j;
	100	int tria_node_ids[3];
	101	int tria_matice_id;
	102	int tria_matpar_id;
	103	double nodeinputs[3];
	104
	105	/id: /
	106	this->id=index+1;
	107
	108	/hooks: /
	109	//go recover node ids, needed to initialize the node hook.
	110	if (iomodel->analysis_type==Prognostic2AnalysisEnum \|\| iomodel->analysis_type==Balancedthickness2AnalysisEnum){
	111	/Discontinuous Galerkin/
	112	tria_node_ids[0]=3*index+1;
	113	tria_node_ids[1]=3*index+2;
	114	tria_node_ids[2]=3*index+3;
	115	}
	116	else{
	117	/Continuous Galerkin/
	118	for(i=0;i<3;i++){
	119	tria_node_ids[i]=(int)(iomodel->elements+3index+i); //ids for vertices are in the elements array from Matlab
	120	}
	121	}
	122	tria_matice_id=index+1; //refers to the corresponding ice material object
	123	tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object
	124
	125	this->hnodes.Init(tria_node_ids,3);
	126	this->hmatice.Init(&tria_matice_id,1);
	127	this->hmatpar.Init(&tria_matpar_id,1);
	128
	129	//intialize inputs, and add as many inputs per element as requested:
	130	this->inputs=new Inputs();
	131
	132	if (iomodel->thickness) {
	133	for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_node_ids[i]-1];
	134	this->inputs->AddInput(new PengridVertexInput(ThicknessEnum,nodeinputs));
	135	}
	136	if (iomodel->surface) {
	137	for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_node_ids[i]-1];
	138	this->inputs->AddInput(new PengridVertexInput(SurfaceEnum,nodeinputs));
	139	}
	140	if (iomodel->bed) {
	141	for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_node_ids[i]-1];
	142	this->inputs->AddInput(new PengridVertexInput(BedEnum,nodeinputs));
	143	}
	144	if (iomodel->drag_coefficient) {
	145	for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_node_ids[i]-1];
	146	this->inputs->AddInput(new PengridVertexInput(DragCoefficientEnum,nodeinputs));
	147
	148	if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index]));
	149	if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index]));
	150	this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type));
	151
	152	}
	153	if (iomodel->melting_rate) {
	154	for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_node_ids[i]-1];
	155	this->inputs->AddInput(new PengridVertexInput(MeltingRateEnum,nodeinputs));
	156	}
	157	if (iomodel->accumulation_rate) {
	158	for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_node_ids[i]-1];
	159	this->inputs->AddInput(new PengridVertexInput(AccumulationRateEnum,nodeinputs));
	160	}
	161	if (iomodel->geothermalflux) {
	162	for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_node_ids[i]-1];
	163	this->inputs->AddInput(new PengridVertexInput(GeothermalFluxEnum,nodeinputs));
	164	}
	165
	166	if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index]));
	167	if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index]));
	168	if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index]));
	169	if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index]));
	170
	171	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
	172	this->parameters=NULL;
	173
	174
	175	}
	176	/}}}/
[2907]	177	/FUNCTION Pengrid::destructor {{{1/
[387]	178	Pengrid::~Pengrid(){
	179	return;
	180	}
[2907]	181	/}}}1/
[387]	182
[2907]	183	/Object marshall/
	184	/FUNCTION Pengrid::Marshall {{{1/
[387]	185	void Pengrid::Marshall(char** pmarshalled_dataset){
	186
	187	char* marshalled_dataset=NULL;
	188	int enum_type=0;
	189
	190	/recover marshalled_dataset: /
	191	marshalled_dataset=*pmarshalled_dataset;
	192
	193	/get enum type of Pengrid: /
[3567]	194	enum_type=PengridEnum;
[387]	195
	196	/marshall enum: /
	197	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
	198
	199	/marshall Pengrid data: /
	200	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
[483]	201	memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
[387]	202	memcpy(marshalled_dataset,&dof,sizeof(dof));marshalled_dataset+=sizeof(dof);
	203	memcpy(marshalled_dataset,&active,sizeof(active));marshalled_dataset+=sizeof(active);
	204	memcpy(marshalled_dataset,&penalty_offset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
	205	memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
	206	memcpy(marshalled_dataset,&node_id,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
	207	memcpy(marshalled_dataset,&node_offset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
[483]	208	memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
	209	memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
[2342]	210	memcpy(marshalled_dataset,&stabilize_constraints,sizeof(stabilize_constraints));marshalled_dataset+=sizeof(stabilize_constraints);
	211	memcpy(marshalled_dataset,&zigzag_counter,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
[387]	212
	213	*pmarshalled_dataset=marshalled_dataset;
	214	return;
	215	}
[2907]	216	/}}}1/
	217	/FUNCTION Pengrid::MarshallSize {{{1/
[387]	218	int Pengrid::MarshallSize(){
	219
	220	return sizeof(id)+
[483]	221	sizeof(mparid)+
[387]	222	sizeof(dof)+
	223	sizeof(active)+
	224	sizeof(penalty_offset)+
	225	sizeof(thermal_steadystate)+
	226	sizeof(node_id)+
	227	sizeof(node_offset)+
[488]	228	sizeof(matpar)+
	229	sizeof(matpar_offset)+
[2342]	230	sizeof(stabilize_constraints)+
	231	sizeof(zigzag_counter)+
[387]	232	sizeof(int); //sizeof(int) for enum type
	233	}
[2907]	234	/}}}1/
	235	/FUNCTION Pengrid::Demarshall {{{1/
[387]	236	void Pengrid::Demarshall(char** pmarshalled_dataset){
	237
	238	char* marshalled_dataset=NULL;
	239
	240	/recover marshalled_dataset: /
	241	marshalled_dataset=*pmarshalled_dataset;
	242
	243	/*this time, no need to get enum type, the pointer directly points to the beginning of the
	244	object data (thanks to DataSet::Demarshall):/
	245
	246	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
[492]	247	memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
[387]	248	memcpy(&dof,marshalled_dataset,sizeof(dof));marshalled_dataset+=sizeof(dof);
	249	memcpy(&active,marshalled_dataset,sizeof(active));marshalled_dataset+=sizeof(active);
	250	memcpy(&penalty_offset,marshalled_dataset,sizeof(penalty_offset));marshalled_dataset+=sizeof(penalty_offset);
	251	memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
	252	memcpy(&node_id,marshalled_dataset,sizeof(node_id));marshalled_dataset+=sizeof(node_id);
	253	memcpy(&node_offset,marshalled_dataset,sizeof(node_offset));marshalled_dataset+=sizeof(node_offset);
[483]	254	memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
	255	memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
[2342]	256	memcpy(&stabilize_constraints,marshalled_dataset,sizeof(stabilize_constraints));marshalled_dataset+=sizeof(stabilize_constraints);
	257	memcpy(&zigzag_counter,marshalled_dataset,sizeof(zigzag_counter));marshalled_dataset+=sizeof(zigzag_counter);
[387]	258
	259	node=NULL;
[483]	260	matpar=NULL;
[387]	261
	262	/return: /
	263	*pmarshalled_dataset=marshalled_dataset;
	264	return;
	265	}
[2907]	266	/}}}1/
[387]	267
[2907]	268	/Object functions/
	269	/FUNCTION Pengrid::copy {{{1/
	270	Object* Pengrid::copy() {
	271	return new Pengrid(*this);
[387]	272	}
[2907]	273	/}}}1/
	274	/FUNCTION Pengrid::Configure {{{1/
[387]	275
	276	void Pengrid::Configure(void* pelementsin,void* pnodesin,void* pmaterialsin){
	277
	278	DataSet* nodesin=NULL;
[483]	279	DataSet* materialsin=NULL;
[387]	280
	281	/Recover pointers :/
	282	nodesin=(DataSet*)pnodesin;
[483]	283	materialsin=(DataSet*)pmaterialsin;
[387]	284
	285	/Link this load with its nodes: /
	286	ResolvePointers((Object**)&node,&node_id,&node_offset,1,nodesin);
[483]	287	ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);
[387]	288	}
[2907]	289	/}}}1/
	290	/FUNCTION Pengrid::CreateKMatrix {{{1/
[387]	291
[3621]	292	void Pengrid::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){
[387]	293
	294	/No loads applied, do nothing: /
	295	return;
	296
	297	}
[2907]	298	/}}}1/
	299	/FUNCTION Pengrid::CreatePVector {{{1/
[3621]	300	void Pengrid::CreatePVector(Vec pg, int analysis_type,int sub_analysis_type){
[387]	301
	302	/No loads applied, do nothing: /
	303	return;
	304
	305	}
[2907]	306	/}}}1/
	307	/FUNCTION Pengrid::DeepEcho {{{1/
	308	void Pengrid::DeepEcho(void){
	309
	310	printf("Pengrid:\n");
	311	printf(" id: %i\n",id);
	312	printf(" mparid: %i\n",mparid);
	313	printf(" dof: %i\n",dof);
	314	printf(" active: %i\n",active);
	315	printf(" penalty_offset: %g\n",penalty_offset);
	316	printf(" thermal_steadystate: %i\n",thermal_steadystate);
	317	printf(" node_id: [%i]\n",node_id);
	318	printf(" node_offset: [%i]\n",node_offset);
	319	printf(" matpar_offset=%i\n",matpar_offset);
[387]	320
[2907]	321	if(node)node->Echo();
	322	if(matpar)matpar->Echo();
	323	return;
	324	}
	325	/}}}1/
	326	/FUNCTION Pengrid::DistributenumDofs {{{1/
	327	void Pengrid::DistributeNumDofs(int* numdofpernode,int analysis_type,int sub_analysis_type){return;}
	328	/}}}1/
	329	/FUNCTION Pengrid::Echo {{{1/
	330	void Pengrid::Echo(void){
	331
	332	printf("Pengrid:\n");
	333	printf(" id: %i\n",id);
	334	printf(" mparid: %i\n",mparid);
	335	printf(" dof: %i\n",dof);
	336	printf(" active: %i\n",active);
	337	printf(" penalty_offset: %g\n",penalty_offset);
	338	printf(" thermal_steadystate: %i\n",thermal_steadystate);
	339	printf(" node_id: [%i]\n",node_id);
	340	printf(" node_offset: [%i]\n",node_offset);
	341	printf(" matpar_offset=%i\n",matpar_offset);
	342
	343	return;
[387]	344	}
[2907]	345	/}}}1/
	346	/FUNCTION Pengrid::Enum {{{1/
	347	int Pengrid::Enum(void){
[387]	348
[3567]	349	return PengridEnum;
[2907]	350	}
	351	/}}}1/
	352	/FUNCTION Pengrid::GetDofList {{{1/
	353	void Pengrid::GetDofList(int* doflist,int* pnumberofdofspernode){
	354
	355	int j;
	356	int doflist_per_node[MAXDOFSPERNODE];
	357	int numberofdofspernode;
	358
	359	node->GetDofList(&doflist_per_node[0],&numberofdofspernode);
	360	for(j=0;j<numberofdofspernode;j++){
	361	doflist[j]=doflist_per_node[j];
	362	}
	363
	364	/Assign output pointers:/
	365	*pnumberofdofspernode=numberofdofspernode;
	366	}
	367	/}}}1/
	368	/FUNCTION Pengrid::GetId {{{1/
	369	int Pengrid::GetId(void){ return id; }
	370	/}}}1/
	371	/FUNCTION Pengrid::GetName {{{1/
	372	char* Pengrid::GetName(void){
	373	return "pengrid";
	374	}
	375	/}}}1/
	376	/FUNCTION Pengrid::MyRank {{{1/
	377	int Pengrid::MyRank(void){
	378	extern int my_rank;
	379	return my_rank;
	380	}
	381	/}}}1/
	382	/FUNCTION Pengrid::PenaltyConstrain {{{1/
[3621]	383	void Pengrid::PenaltyConstrain(int* punstable,int analysis_type,int sub_analysis_type){
[2907]	384
[3567]	385	if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){
[2907]	386
	387	/No penalty to check/
	388	return;
	389
	390	}
[3567]	391	else if (analysis_type==ThermalAnalysisEnum){
[2907]	392
[3621]	393	PenaltyConstrainThermal(punstable,analysis_type,sub_analysis_type);
[2907]	394
	395	}
[3567]	396	else if (analysis_type==MeltingAnalysisEnum){
[2907]	397
	398	/No penalty to check/
	399	return;
	400
	401	}
	402	else{
[3570]	403	ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
[2907]	404	}
	405
	406	}
	407	/}}}1/
	408	/FUNCTION Pengrid::PenaltyConstrainThermal {{{1/
[3621]	409	void Pengrid::PenaltyConstrainThermal(int* punstable,int analysis_type,int sub_analysis_type){
[2907]	410
	411	// The penalty is stable if it doesn't change during to successive iterations.
	412
	413	int found=0;
	414	const int numgrids=1;
	415
	416
	417	double pressure;
	418	double temperature;
	419	double beta,t_pmp;
	420	double meltingpoint;
	421	int new_active;
	422	int dofs1[1]={0};
	423	int unstable=0;
	424	int reset_penalties=0;
	425
	426	ParameterInputs* inputs=NULL;
	427
	428	/check that pengrid is not a clone (penalty to be added only once)/
	429	if (node->IsClone()){
	430	unstable=0;
	431	*punstable=unstable;
	432	return;
	433	}
	434
	435	/recover pointers: /
	436	inputs=(ParameterInputs*)vinputs;
	437
	438	//First recover beta, pressure and temperature vectors;
	439	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
[3332]	440	if(!found)ISSMERROR(" could not find pressure in inputs!");
[2907]	441
	442	found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
[3332]	443	if(!found)ISSMERROR(" could not find temperature in inputs!");
[2907]	444
	445	found=inputs->Recover("reset_penalties",&reset_penalties);
	446
	447	if(reset_penalties)zigzag_counter=0;
	448
	449	//Compute pressure melting point
	450	meltingpoint=matpar->GetMeltingPoint();
	451	beta=matpar->GetBeta();
	452
	453	t_pmp=meltingpoint-beta*pressure;
	454
	455	//Figure out if temperature is over melting_point, in which case, this penalty needs to be activated.
	456
	457	if (temperature>t_pmp){
	458	new_active=1;
	459	}
	460	else{
	461	new_active=0;
	462	}
	463
	464
	465	//Figure out stability of this penalty
	466	if (active==new_active){
	467	unstable=0;
	468	}
	469	else{
	470	unstable=1;
	471	if(stabilize_constraints)zigzag_counter++;
	472	}
	473
	474	/If penalty keeps zigzagging more than 5 times: /
	475	if(stabilize_constraints){
	476	if(zigzag_counter>stabilize_constraints){
	477	unstable=0;
	478	active=1;
	479	}
	480	}
	481
	482	//Set penalty flag
	483	active=new_active;
	484
	485	//Assign output pointers:/
	486	*punstable=unstable;
	487	}
	488	/}}}1/
	489	/FUNCTION Pengrid::PenaltyCreateMatrix {{{1/
[3621]	490	void Pengrid::PenaltyCreateKMatrix(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
[394]	491
[3567]	492	if ((analysis_type==DiagnosticAnalysisEnum) && ((sub_analysis_type==StokesAnalysisEnum))){
[394]	493
[3621]	494	PenaltyCreateKMatrixDiagnosticStokes( Kgg,kmax,analysis_type,sub_analysis_type);
[483]	495	}
[3567]	496	else if (analysis_type==ThermalAnalysisEnum){
[483]	497
[3621]	498	PenaltyCreateKMatrixThermal( Kgg,kmax,analysis_type,sub_analysis_type);
[483]	499
	500	}
[3567]	501	else if (analysis_type==MeltingAnalysisEnum){
[483]	502
[3621]	503	PenaltyCreateKMatrixMelting( Kgg,kmax,analysis_type,sub_analysis_type);
[394]	504
	505	}
	506	else{
[3570]	507	ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
[394]	508	}
	509
[387]	510	}
[2907]	511	/}}}1/
	512	/FUNCTION Pengrid::PenaltyCreateKMatrixDiagnosticStokes {{{1/
[3621]	513	void Pengrid::PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
[394]	514
	515	const int numgrids=1;
[461]	516	const int NDOF4=4;
	517	const int numdof=numgrids*NDOF4;
[394]	518	int doflist[numdof];
	519	int numberofdofspernode;
	520
[483]	521	int dofs1[1]={0};
	522	int dofs2[1]={1};
[394]	523	double slope[2];
	524	int found=0;
[442]	525	double Ke[4][4]={0.0};
[394]	526
	527	ParameterInputs* inputs=NULL;
	528
	529	/recover pointers: /
	530	inputs=(ParameterInputs*)vinputs;
	531
	532	/Get dof list: /
	533	GetDofList(&doflist[0],&numberofdofspernode);
	534
	535	/recover slope: /
[483]	536	found=inputs->Recover("bedslopex",&slope[0],1,dofs1,numgrids,(void**)&node);
[3332]	537	if(!found)ISSMERROR(" bedslopex needed in inputs!");
[483]	538	found=inputs->Recover("bedslopey",&slope[1],1,dofs2,numgrids,(void**)&node);
[3332]	539	if(!found)ISSMERROR(" bedslopey needed in inputs!");
[394]	540
	541	//Create elementary matrix: add penalty to contrain wb (wb=ubdb/dx+vbdb/dy)
[1439]	542	Ke[2][0]=-slope[0]kmaxpow((double)10.0,penalty_offset);
	543	Ke[2][1]=-slope[1]kmaxpow((double)10.0,penalty_offset);
	544	Ke[2][2]=kmax*pow((double)10,penalty_offset);
[394]	545
	546	/Add Ke to global matrix Kgg: /
	547	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
	548	}
[2907]	549	/}}}1/
	550	/FUNCTION Pengrid::PenaltyCreateKMatrixMelting {{{1/
[3621]	551	void Pengrid::PenaltyCreateKMatrixMelting(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
[483]	552
[926]	553
[483]	554	int found=0;
	555	const int numgrids=1;
	556	const int NDOF1=1;
	557	const int numdof=numgrids*NDOF1;
	558	double Ke[numdof][numdof]={0.0};
	559	int dofs1[1]={0};
	560	int doflist[numdof];
	561	int numberofdofspernode;
	562	double meltingpoint;
	563
	564	double pressure;
	565	double temperature;
	566	double beta,t_pmp;
	567
	568	ParameterInputs* inputs=NULL;
	569
[926]	570	/check that pengrid is not a clone (penalty to be added only once)/
	571	if (node->IsClone()) return;
	572
[483]	573	/recover pointers: /
	574	inputs=(ParameterInputs*)vinputs;
	575
	576	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
[3332]	577	if(!found)ISSMERROR(" could not find pressure in inputs!");
[483]	578
	579	found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
[3332]	580	if(!found)ISSMERROR(" could not find temperature in inputs!");
[483]	581
	582	/Get dof list: /
	583	GetDofList(&doflist[0],&numberofdofspernode);
	584
	585	//Compute pressure melting point
	586	meltingpoint=matpar->GetMeltingPoint();
	587	beta=matpar->GetBeta();
	588	t_pmp=meltingpoint-beta*pressure;
	589
	590	//Add penalty load
	591	if (temperature<t_pmp){ //If T<Tpmp, there must be no melting. Therefore, melting should be constrained to 0 when T<Tpmp, instead of using spcs, use penalties
[1439]	592	Ke[0][0]=kmax*pow((double)10,penalty_offset);
[483]	593	}
	594
	595	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
	596	}
[2907]	597	/}}}1/
	598	/FUNCTION Pengrid::PenaltyCreateKMatrixThermal {{{1/
[3621]	599	void Pengrid::PenaltyCreateKMatrixThermal(Mat Kgg,double kmax,int analysis_type,int sub_analysis_type){
[483]	600
[2907]	601	int found=0;
	602
	603	const int numgrids=1;
	604	const int NDOF1=1;
	605	const int numdof=numgrids*NDOF1;
	606	double Ke[numdof][numdof];
	607	int doflist[numdof];
	608	int numberofdofspernode;
	609
	610	ParameterInputs* inputs=NULL;
	611
	612	/recover pointers: /
	613	inputs=(ParameterInputs*)vinputs;
	614
	615
	616	if(!active)return;
	617
	618	/Get dof list: /
	619	GetDofList(&doflist[0],&numberofdofspernode);
	620
	621	Ke[0][0]=kmax*pow((double)10,penalty_offset);
	622
	623	/Add Ke to global matrix Kgg: /
	624	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke,ADD_VALUES);
	625	}
	626	/}}}1/
	627	/FUNCTION Pengrid::PenaltyCreatePVector {{{1/
[3621]	628	void Pengrid::PenaltyCreatePVector(Vec pg,double kmax,int analysis_type,int sub_analysis_type){
[442]	629
[3567]	630	if (analysis_type==ThermalAnalysisEnum){
[483]	631
[3621]	632	PenaltyCreatePVectorThermal( pg,kmax,analysis_type,sub_analysis_type);
[483]	633
	634	}
[3567]	635	else if (analysis_type==MeltingAnalysisEnum){
[483]	636
[3621]	637	PenaltyCreatePVectorMelting( pg,kmax,analysis_type,sub_analysis_type);
[442]	638
[483]	639	}
[3567]	640	else if (analysis_type==DiagnosticAnalysisEnum){
[498]	641
	642	/No loads applied, do nothing: /
	643	return;
	644
	645	}
[483]	646	else{
[3570]	647	ISSMERROR("%s%i%s%i%s","analysis: ",analysis_type," and sub_analysis_type: ",sub_analysis_type," not supported yet");
[483]	648	}
	649
[387]	650	}
[2907]	651	/}}}1/
	652	/FUNCTION Pengrid::PenaltyCreatePVectorMelting {{{1/
[3621]	653	void Pengrid::PenaltyCreatePVectorMelting(Vec pg, double kmax,int analysis_type,int sub_analysis_type){
[483]	654
	655	const int numgrids=1;
	656	const int NDOF1=1;
	657	const int numdof=numgrids*NDOF1;
[539]	658	int doflist[numdof];
	659	double P_terms[numdof]={0.0};
	660	int numberofdofspernode;
[483]	661	int found=0;
[539]	662	int dofs1[1]={0};
[483]	663	double pressure;
	664	double temperature;
	665	double melting_offset;
	666	double meltingpoint;
	667	double beta, heatcapacity;
	668	double latentheat;
	669	double t_pmp;
	670	double dt;
	671
	672	ParameterInputs* inputs=NULL;
	673
[926]	674	/check that pengrid is not a clone (penalty to be added only once)/
	675	if (node->IsClone()) return;
	676
[483]	677	/recover pointers: /
	678	inputs=(ParameterInputs*)vinputs;
	679
[539]	680	/Get dof list: /
	681	GetDofList(&doflist[0],&numberofdofspernode);
	682
[483]	683	//First recover pressure,melting offset and temperature vectors
	684	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
[3332]	685	if(!found)ISSMERROR(" could not find pressure in inputs!");
[483]	686
	687	found=inputs->Recover("temperature",&temperature,1,dofs1,numgrids,(void**)&node);
[3332]	688	if(!found)ISSMERROR(" could not find temperature in inputs!");
[483]	689
	690	found=inputs->Recover("melting_offset",&melting_offset);
[3332]	691	if(!found)ISSMERROR(" could not find melting_offset in inputs!");
[483]	692
	693	found=inputs->Recover("dt",&dt);
[3332]	694	if(!found)ISSMERROR(" could not find dt in inputs!");
[483]	695
	696	meltingpoint=matpar->GetMeltingPoint();
	697	beta=matpar->GetBeta();
	698	heatcapacity=matpar->GetHeatCapacity();
	699	latentheat=matpar->GetLatentHeat();
	700
	701	//Compute pressure melting point
	702	t_pmp=meltingpoint-beta*pressure;
	703
	704	//Add penalty load
	705	//This time, the penalty must have the same value as the one used for the thermal computation
	706	//so that the corresponding melting can be computed correctly
	707	//In the thermal computation, we used kmax=melting_offset, and the same penalty_offset
	708	if (temperature<t_pmp){ //%no melting
	709	P_terms[0]=0;
	710	}
	711	else{
[1901]	712	if (dt){
	713	P_terms[0]=melting_offsetpow((double)10,penalty_offset)(temperature-t_pmp)/dt;
[483]	714	}
	715	else{
[1901]	716	P_terms[0]=melting_offsetpow((double)10,penalty_offset)(temperature-t_pmp);
[483]	717	}
	718	}
[539]	719
[483]	720	/Add P_terms to global vector pg: /
	721	VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
	722	}
[2907]	723	/}}}1/
	724	/FUNCTION Pengrid::PenaltyCreatePVectorThermal {{{1/
[3621]	725	void Pengrid::PenaltyCreatePVectorThermal(Vec pg, double kmax,int analysis_type,int sub_analysis_type){
[483]	726
	727	const int numgrids=1;
[2907]	728	const int NDOF1=1;
	729	const int numdof=numgrids*NDOF1;
	730	int doflist[numdof];
	731	double P_terms[numdof]={0.0};
	732	int numberofdofspernode;
	733	int found=0;
[483]	734	double pressure;
[2907]	735	int dofs1[1]={0};
[483]	736	double meltingpoint;
[2907]	737	double beta;
	738	double t_pmp;
	739
[483]	740	ParameterInputs* inputs=NULL;
	741
	742	/recover pointers: /
	743	inputs=(ParameterInputs*)vinputs;
	744
[2907]	745	if(!active)return;
	746
	747	/Get dof list: /
	748	GetDofList(&doflist[0],&numberofdofspernode);
	749
	750	//First recover pressure
[483]	751	found=inputs->Recover("pressure",&pressure,1,dofs1,numgrids,(void**)&node);
[3332]	752	if(!found)ISSMERROR(" could not find pressure in inputs!");
[483]	753
	754	//Compute pressure melting point
	755	meltingpoint=matpar->GetMeltingPoint();
	756	beta=matpar->GetBeta();
	757	t_pmp=meltingpoint-beta*pressure;
	758
[2907]	759	//Add penalty load
	760	P_terms[0]=kmaxpow((double)10,penalty_offset)t_pmp;
[483]	761
[2907]	762	/Add P_terms to global vector pg: /
	763	VecSetValues(pg,numdof,doflist,(const double*)P_terms,ADD_VALUES);
[483]	764	}
[2907]	765	/}}}1/
	766	/FUNCTION Pengrid::UpdateFromInputs {{{1/
	767	void Pengrid::UpdateFromInputs(void* inputs){
	768
	769	}
	770	/}}}1/

Note: See TracBrowser for help on using the repository browser.

Download in other formats: