source: issm/oecreview/Archive/18296-19100/ISSM-18343-18344.diff

Last change on this file was 19102, checked in by Mathieu Morlighem, 10 years ago

NEW: added 18296-19100
File size: 24.0 KB
RevLine 
[19102]1Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
2===================================================================
3--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18343)
4+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 18344)
5@@ -70,6 +70,9 @@
6 if(stabilization==3){
7 iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for DG, we need the spc in the element
8 }
9+ if(stabilization==4){
10+ iomodel->FetchDataToInput(elements,MasstransportSpcthicknessEnum); //for FCT, we need the spc in the element (penlaties)
11+ }
12
13 if(iomodel->domaintype!=Domain2DhorizontalEnum){
14 iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
15@@ -147,10 +150,11 @@
16 int stabilization;
17 iomodel->Constant(&stabilization,MasstransportStabilizationEnum);
18
19- /*Do not add constraints in DG, they are weakly imposed*/
20+ /*Do not add constraints in DG, they are weakly imposed*/
21 if(stabilization!=3){
22 IoModelToConstraintsx(constraints,iomodel,MasstransportSpcthicknessEnum,MasstransportAnalysisEnum,P1Enum);
23 }
24+ /*Do not add constraints in FCT, they are imposed using penalties*/
25 }/*}}}*/
26 void MasstransportAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
27
28@@ -674,7 +678,7 @@
29
30 }/*}}}*/
31 void MasstransportAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
32- _error_("not implemented yet");
33+ element->GetSolutionFromInputsOneDof(solution,ThicknessEnum);
34 }/*}}}*/
35 void MasstransportAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
36 _error_("Not implemented yet");
37@@ -772,3 +776,148 @@
38 }
39 return;
40 }/*}}}*/
41+
42+/*Flux Correction Transport*/
43+void MasstransportAnalysis::LumpedMassMatrix(Vector<IssmDouble>** pMlff,FemModel* femmodel){/*{{{*/
44+
45+ /*Intermediaries*/
46+ int configuration_type;
47+
48+ /*Initialize Lumped mass matrix (actually we just save its diagonal)*/
49+ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
50+ int fsize = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum);
51+ int flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum);
52+ Vector<IssmDouble>* Mlff = new Vector<IssmDouble>(flocalsize,fsize);
53+
54+ /*Create and assemble matrix*/
55+ for(int i=0;i<femmodel->elements->Size();i++){
56+ Element* element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
57+ ElementMatrix* MLe = this->CreateMassMatrix(element);
58+ if(MLe){
59+ MLe->Lump();
60+ MLe->AddDiagonalToGlobal(Mlff);
61+ }
62+ delete MLe;
63+ }
64+ Mlff->Assemble();
65+
66+ /*Assign output pointer*/
67+ *pMlff=Mlff;
68+}/*}}}*/
69+ElementMatrix* MasstransportAnalysis::CreateMassMatrix(Element* element){/*{{{*/
70+
71+ /* Check if ice in element */
72+ if(!element->IsIceInElement()) return NULL;
73+
74+ /*Intermediaries*/
75+ IssmDouble D,Jdet;
76+ IssmDouble* xyz_list = NULL;
77+
78+ /*Fetch number of nodes and dof for this finite element*/
79+ int numnodes = element->GetNumberOfNodes();
80+
81+ /*Initialize Element vector and other vectors*/
82+ ElementMatrix* Me = element->NewElementMatrix();
83+ IssmDouble* basis = xNew<IssmDouble>(numnodes);
84+
85+ /*Retrieve all inputs and parameters*/
86+ element->GetVerticesCoordinates(&xyz_list);
87+
88+ /* Start looping on the number of gaussian points: */
89+ Gauss* gauss=element->NewGauss(2);
90+ for(int ig=gauss->begin();ig<gauss->end();ig++){
91+ gauss->GaussPoint(ig);
92+
93+ element->JacobianDeterminant(&Jdet,xyz_list,gauss);
94+ element->NodalFunctions(basis,gauss);
95+
96+ D=gauss->weight*Jdet;
97+ TripleMultiply(basis,1,numnodes,1,
98+ &D,1,1,0,
99+ basis,1,numnodes,0,
100+ &Me->values[0],1);
101+ }
102+
103+ /*Clean up and return*/
104+ xDelete<IssmDouble>(xyz_list);
105+ xDelete<IssmDouble>(basis);
106+ delete gauss;
107+ return Me;
108+}/*}}}*/
109+void MasstransportAnalysis::FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel){/*{{{*/
110+
111+ /*Output*/
112+ Matrix<IssmDouble>* Kff = NULL;
113+ Matrix<IssmDouble>* Kfs = NULL;
114+
115+ /*Initialize Jacobian Matrix*/
116+ AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel);
117+
118+ /*Create and assemble matrix*/
119+ for(int i=0;i<femmodel->elements->Size();i++){
120+ Element* element = dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
121+ ElementMatrix* Ke = this->CreateFctKMatrix(element);
122+ if(Ke) Ke->AddToGlobal(Kff,Kfs);
123+ delete Ke;
124+ }
125+ Kff->Assemble();
126+ Kfs->Assemble();
127+
128+ /*Assign output pointer*/
129+ *pKff=Kff;
130+ *pKfs=Kfs;
131+}/*}}}*/
132+ElementMatrix* MasstransportAnalysis::CreateFctKMatrix(Element* element){/*{{{*/
133+
134+ /* Check if ice in element */
135+ if(!element->IsIceInElement()) return NULL;
136+
137+ /*Intermediaries */
138+ IssmDouble Jdet;
139+ IssmDouble vx,vy;
140+ IssmDouble* xyz_list = NULL;
141+
142+ /*Fetch number of nodes and dof for this finite element*/
143+ int numnodes = element->GetNumberOfNodes();
144+ int dim = 2;
145+
146+ /*Initialize Element vector and other vectors*/
147+ ElementMatrix* Ke = element->NewElementMatrix();
148+ IssmDouble* B = xNew<IssmDouble>(dim*numnodes);
149+ IssmDouble* Bprime = xNew<IssmDouble>(dim*numnodes);
150+ IssmDouble* D = xNewZeroInit<IssmDouble>(dim*dim);
151+
152+ /*Retrieve all inputs and parameters*/
153+ element->GetVerticesCoordinates(&xyz_list);
154+ Input* vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input);
155+ Input* vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input);
156+
157+ /* Start looping on the number of gaussian points: */
158+ Gauss* gauss=element->NewGauss(2);
159+ for(int ig=gauss->begin();ig<gauss->end();ig++){
160+ gauss->GaussPoint(ig);
161+
162+ element->JacobianDeterminant(&Jdet,xyz_list,gauss);
163+ GetB(B,element,dim,xyz_list,gauss);
164+ GetBprime(Bprime,element,dim,xyz_list,gauss);
165+ vxaverage_input->GetInputValue(&vx,gauss);
166+ vyaverage_input->GetInputValue(&vy,gauss);
167+
168+ D[0*dim+0] = -gauss->weight*vx*Jdet;
169+ D[1*dim+1] = -gauss->weight*vy*Jdet;
170+
171+ TripleMultiply(B,dim,numnodes,1,
172+ D,dim,dim,0,
173+ Bprime,dim,numnodes,0,
174+ &Ke->values[0],1);
175+
176+ }
177+
178+ /*Clean up and return*/
179+ xDelete<IssmDouble>(xyz_list);
180+ xDelete<IssmDouble>(B);
181+ xDelete<IssmDouble>(Bprime);
182+ xDelete<IssmDouble>(D);
183+ delete gauss;
184+ return Ke;
185+}/*}}}*/
186Index: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.h
187===================================================================
188--- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 18343)
189+++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.h (revision 18344)
190@@ -35,5 +35,11 @@
191 void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
192 void InputUpdateFromSolution(IssmDouble* solution,Element* element);
193 void UpdateConstraints(FemModel* femmodel);
194+
195+ /*FCT*/
196+ void LumpedMassMatrix(Vector<IssmDouble>** pMLff,FemModel* femmodel);
197+ void FctKMatrix(Matrix<IssmDouble>** pKff,Matrix<IssmDouble>** pKfs,FemModel* femmodel);
198+ ElementMatrix* CreateMassMatrix(Element* element);
199+ ElementMatrix* CreateFctKMatrix(Element* element);
200 };
201 #endif
202Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp
203===================================================================
204--- ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 0)
205+++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 18344)
206@@ -0,0 +1,161 @@
207+/*!\file: solutionsequence_linear.cpp
208+ * \brief: numerical core of linear solutions
209+ */
210+
211+#include "../toolkits/toolkits.h"
212+#include "../classes/classes.h"
213+#include "../shared/shared.h"
214+#include "../modules/modules.h"
215+
216+void solutionsequence_fct(FemModel* femmodel){
217+
218+ /*intermediary: */
219+ Vector<IssmDouble>* Mlf = NULL;
220+ Matrix<IssmDouble>* Kff = NULL;
221+ Matrix<IssmDouble>* Kfs = NULL;
222+ Vector<IssmDouble>* ug = NULL;
223+ Vector<IssmDouble>* uf = NULL;
224+ Vector<IssmDouble>* ys = NULL;
225+
226+ IssmDouble theta,deltat;
227+ int dof,ncols,ncols2,ncols3,rstart,rend;
228+ int configuration_type,analysis_type;
229+ double d,diagD,mi;
230+ double dmax = 0.;
231+ int* cols = NULL;
232+ int* cols2 = NULL;
233+ int* cols3 = NULL;
234+ double* vals = NULL;
235+ double* vals2 = NULL;
236+ double* vals3 = NULL;
237+
238+ /*Create analysis*/
239+ MasstransportAnalysis* analysis = new MasstransportAnalysis();
240+
241+ /*Recover parameters: */
242+ femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
243+ femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
244+ femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
245+ femmodel->UpdateConstraintsx();
246+ theta = 0.5;
247+
248+ /*Create lumped mass matrix*/
249+ analysis->LumpedMassMatrix(&Mlf,femmodel);
250+ analysis->FctKMatrix(&Kff,&Kfs,femmodel);
251+
252+ /*Create Dff Matrix*/
253+ #ifdef _HAVE_PETSC_
254+ Mat Kff_transp = NULL;
255+ Mat Dff_petsc = NULL;
256+ Mat LHS = NULL;
257+ Mat Kff_petsc = Kff->pmatrix->matrix;
258+ Vec Mlf_petsc = Mlf->pvector->vector;
259+ MatTranspose(Kff_petsc,MAT_INITIAL_MATRIX,&Kff_transp);
260+ MatDuplicate(Kff_petsc,MAT_SHARE_NONZERO_PATTERN,&Dff_petsc);
261+ MatGetOwnershipRange(Kff_transp,&rstart,&rend);
262+ for(int row=rstart; row<rend; row++){
263+ diagD = 0.;
264+ MatGetRow(Kff_petsc ,row,&ncols, (const int**)&cols, (const double**)&vals);
265+ MatGetRow(Kff_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
266+ _assert_(ncols==ncols2);
267+ for(int j=0; j<ncols; j++) {
268+ _assert_(cols[j]==cols2[j]);
269+ d = max(max(-vals[j],-vals2[j]),0.);
270+ MatSetValues(Dff_petsc,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
271+ if(cols[j]!=row) diagD -= d;
272+ }
273+ MatSetValues(Dff_petsc,1,&row,1,&row,(const double*)&diagD,INSERT_VALUES);
274+ MatRestoreRow(Kff_petsc, row,&ncols, (const int**)&cols, (const double**)&vals);
275+ MatRestoreRow(Kff_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
276+ }
277+ MatAssemblyBegin(Dff_petsc,MAT_FINAL_ASSEMBLY);
278+ MatAssemblyEnd( Dff_petsc,MAT_FINAL_ASSEMBLY);
279+ MatFree(&Kff_transp);
280+
281+ /*Create LHS: [ML − theta*detlat *(K+D)^n+1]*/
282+ MatDuplicate(Kff_petsc,MAT_SHARE_NONZERO_PATTERN,&LHS);
283+ for(int row=rstart; row<rend; row++){
284+ MatGetRow(Kff_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
285+ MatGetRow(Dff_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
286+ _assert_(ncols==ncols2);
287+ for(int j=0; j<ncols; j++) {
288+ _assert_(cols[j]==cols2[j]);
289+ d = -theta*deltat*(vals[j] + vals2[j]);
290+ if(cols[j]==row){
291+ VecGetValues(Mlf_petsc,1,(const int*)&cols[j],&mi);
292+ d += mi;
293+ }
294+ if(fabs(d)>dmax) dmax = fabs(d);
295+ MatSetValues(LHS,1,&row,1,&cols[j],(const double*)&d,INSERT_VALUES);
296+ }
297+ MatRestoreRow(Kff_petsc,row,&ncols, (const int**)&cols, (const double**)&vals);
298+ MatRestoreRow(Dff_petsc,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
299+ }
300+
301+ /*Penalize Dirichlet boundary*/
302+ dmax = dmax * 1.e+3;
303+ for(int i=0;i<femmodel->constraints->Size();i++){
304+ Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
305+ if(constraint->InAnalysis(analysis_type)){
306+ constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
307+ if(dof!=-1){
308+ MatSetValues(LHS,1,&dof,1,&dof,(const double*)&dmax,INSERT_VALUES);
309+ }
310+ }
311+ }
312+ MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
313+ MatAssemblyEnd( LHS,MAT_FINAL_ASSEMBLY);
314+
315+ /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
316+ GetSolutionFromInputsx(&ug,femmodel);
317+ Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
318+ delete ug;
319+ Vec u = uf->pvector->vector;
320+ Vec Ku = NULL;
321+ Vec Du = NULL;
322+ Vec RHS = NULL;
323+ VecDuplicate(u,&Ku);
324+ VecDuplicate(u,&Du);
325+ VecDuplicate(u,&RHS);
326+ MatMult(Kff_petsc,u,Ku);
327+ MatMult(Dff_petsc,u,Du);
328+ VecPointwiseMult(RHS,Mlf_petsc,u);
329+ VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);// RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u
330+ VecFree(&Ku);
331+ VecFree(&Du);
332+ MatFree(&Dff_petsc);
333+ delete uf;
334+
335+ /*Penalize Dirichlet boundary*/
336+ for(int i=0;i<femmodel->constraints->Size();i++){
337+ Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
338+ if(constraint->InAnalysis(analysis_type)){
339+ constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
340+ d = d*dmax;
341+ if(dof!=-1){
342+ VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
343+ }
344+ }
345+ }
346+
347+ /*Go solve!*/
348+ SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
349+ MatFree(&LHS);
350+ VecFree(&RHS);
351+ uf =new Vector<IssmDouble>(u);
352+ VecFree(&u);
353+
354+ Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
355+ InputUpdateFromSolutionx(femmodel,ug);
356+ delete ug;
357+
358+ #else
359+ _error_("PETSc needs to be installed");
360+ #endif
361+
362+ delete Mlf;
363+ delete Kff;
364+ delete Kfs;
365+ delete analysis;
366+
367+}
368Index: ../trunk-jpl/src/c/solutionsequences/solutionsequences.h
369===================================================================
370--- ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 18343)
371+++ ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 18344)
372@@ -15,6 +15,7 @@
373 void solutionsequence_hydro_nonlinear(FemModel* femmodel);
374 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
375 void solutionsequence_newton(FemModel* femmodel);
376+void solutionsequence_fct(FemModel* femmodel);
377 void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
378 void solutionsequence_linear(FemModel* femmodel);
379 void solutionsequence_la(FemModel* femmodel);
380Index: ../trunk-jpl/src/c/Makefile.am
381===================================================================
382--- ../trunk-jpl/src/c/Makefile.am (revision 18343)
383+++ ../trunk-jpl/src/c/Makefile.am (revision 18344)
384@@ -353,6 +353,7 @@
385 ./solutionsequences/solutionsequence_linear.cpp\
386 ./solutionsequences/solutionsequence_nonlinear.cpp\
387 ./solutionsequences/solutionsequence_newton.cpp\
388+ ./solutionsequences/solutionsequence_fct.cpp\
389 ./solutionsequences/convergence.cpp\
390 ./classes/Options/Options.h\
391 ./classes/Options/Options.cpp\
392Index: ../trunk-jpl/src/c/cores/masstransport_core.cpp
393===================================================================
394--- ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 18343)
395+++ ../trunk-jpl/src/c/cores/masstransport_core.cpp (revision 18344)
396@@ -15,7 +15,7 @@
397 int numoutputs,domaintype;
398 bool save_results;
399 bool isFS,isfreesurface,dakota_analysis;
400- int solution_type;
401+ int solution_type,stabilization;
402 char** requested_outputs = NULL;
403
404 /*activate configuration*/
405@@ -28,6 +28,7 @@
406 femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
407 femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
408 femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
409+ femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
410 if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
411
412 /*Calculate new Surface Mass Balance (SMB)*/
413@@ -48,7 +49,12 @@
414 }
415 else{
416 if(VerboseSolution()) _printf0_(" call computational core\n");
417- solutionsequence_linear(femmodel);
418+ if(stabilization==4){
419+ solutionsequence_fct(femmodel);
420+ }
421+ else{
422+ solutionsequence_linear(femmodel);
423+ }
424 femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
425 extrudefrombase_core(femmodel);
426 femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
427Index: ../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp
428===================================================================
429--- ../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp (revision 18343)
430+++ ../trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp (revision 18344)
431@@ -328,6 +328,39 @@
432
433 }
434 /*}}}*/
435+void ElementMatrix::AddDiagonalToGlobal(Vector<IssmDouble>* pf){/*{{{*/
436+
437+ IssmDouble* localvalues=NULL;
438+
439+ /*Check that pf is not NULL*/
440+ _assert_(pf);
441+
442+ /*In debugging mode, check consistency (no NaN, and values not too big)*/
443+ this->CheckConsistency();
444+
445+ if(this->dofsymmetrical){
446+ /*only use row dofs to add values into global matrices: */
447+
448+ if(this->row_fsize){
449+ /*first, retrieve values that are in the f-set from the g-set values matrix: */
450+ localvalues=xNew<IssmDouble>(this->row_fsize);
451+ for(int i=0;i<this->row_fsize;i++){
452+ localvalues[i] = this->values[this->ncols*this->row_flocaldoflist[i]+ this->row_flocaldoflist[i]];
453+ }
454+
455+ /*add local values into global matrix, using the fglobaldoflist: */
456+ pf->SetValues(this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
457+
458+ /*Free ressources:*/
459+ xDelete<IssmDouble>(localvalues);
460+ }
461+ }
462+ else{
463+ _error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
464+ }
465+
466+}
467+/*}}}*/
468 void ElementMatrix::CheckConsistency(void){/*{{{*/
469 /*Check element matrix values, only in debugging mode*/
470 #ifdef _ISSM_DEBUG_
471Index: ../trunk-jpl/src/c/classes/matrix/ElementMatrix.h
472===================================================================
473--- ../trunk-jpl/src/c/classes/matrix/ElementMatrix.h (revision 18343)
474+++ ../trunk-jpl/src/c/classes/matrix/ElementMatrix.h (revision 18344)
475@@ -59,6 +59,7 @@
476 /*ElementMatrix specific routines*/
477 void AddToGlobal(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
478 void AddToGlobal(Matrix<IssmDouble>* Jff);
479+ void AddDiagonalToGlobal(Vector<IssmDouble>* pf);
480 void Echo(void);
481 void CheckConsistency(void);
482 void StaticCondensation(int numindices,int* indices);
483Index: ../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h
484===================================================================
485--- ../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h (revision 18343)
486+++ ../trunk-jpl/src/c/classes/Constraints/SpcDynamic.h (revision 18344)
487@@ -37,6 +37,7 @@
488 /*Constraint virtual functions definitions: {{{*/
489 void ConstrainNode(Nodes* nodes,Parameters* parameters);
490 bool InAnalysis(int analysis_type);
491+ void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
492 /*}}}*/
493 /*SpcDynamic management:{{{ */
494 int GetNodeId();
495Index: ../trunk-jpl/src/c/classes/Constraints/SpcStatic.h
496===================================================================
497--- ../trunk-jpl/src/c/classes/Constraints/SpcStatic.h (revision 18343)
498+++ ../trunk-jpl/src/c/classes/Constraints/SpcStatic.h (revision 18344)
499@@ -36,6 +36,7 @@
500 /*Constraint virtual functions definitions: {{{*/
501 void ConstrainNode(Nodes* nodes,Parameters* parameters);
502 bool InAnalysis(int analysis_type);
503+ void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters){_error_("not implemented yet");};
504 /*}}}*/
505 /*SpcStatic management:{{{ */
506 int GetNodeId();
507Index: ../trunk-jpl/src/c/classes/Constraints/SpcTransient.h
508===================================================================
509--- ../trunk-jpl/src/c/classes/Constraints/SpcTransient.h (revision 18343)
510+++ ../trunk-jpl/src/c/classes/Constraints/SpcTransient.h (revision 18344)
511@@ -13,13 +13,14 @@
512 class SpcTransient: public Constraint{
513
514 private:
515- int sid; /*! id, to track it*/
516- int nodeid; /*!node id*/
517- int dof; /*!component*/
518- IssmDouble* values; /*different values in time*/
519- IssmDouble* times; /*different time steps*/
520- int nsteps; /*number of time steps*/
521- int analysis_type;
522+ bool penalty; /*Is this a penalty constraint */
523+ int sid; /* id, to track it */
524+ int nodeid; /*node id */
525+ int dof; /*component */
526+ IssmDouble *values; /*different values in time */
527+ IssmDouble *times; /*different time steps */
528+ int nsteps; /*number of time steps */
529+ int analysis_type;
530
531 public:
532
533@@ -38,6 +39,7 @@
534 /*Constraint virtual functions definitions: {{{*/
535 void ConstrainNode(Nodes* nodes,Parameters* parameters);
536 bool InAnalysis(int analysis_type);
537+ void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters);
538 /*}}}*/
539 /*SpcTransient management:{{{ */
540 int GetNodeId();
541Index: ../trunk-jpl/src/c/classes/Constraints/Constraint.h
542===================================================================
543--- ../trunk-jpl/src/c/classes/Constraints/Constraint.h (revision 18343)
544+++ ../trunk-jpl/src/c/classes/Constraints/Constraint.h (revision 18344)
545@@ -20,6 +20,7 @@
546
547 virtual ~Constraint(){};
548 virtual void ConstrainNode(Nodes* nodes,Parameters* parameters)=0;
549+ virtual void PenaltyDofAndValue(int* dof,IssmDouble* value,Nodes* nodes,Parameters* parameters)=0;
550 virtual bool InAnalysis(int analysis_type)=0;
551
552 };
553Index: ../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp
554===================================================================
555--- ../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp (revision 18343)
556+++ ../trunk-jpl/src/c/classes/Constraints/SpcTransient.cpp (revision 18344)
557@@ -14,25 +14,27 @@
558
559 /*SpcTransient constructors and destructor*/
560 SpcTransient::SpcTransient(){/*{{{*/
561- sid=-1;
562- nodeid=-1;
563- dof=-1;
564- values=NULL;
565- times=NULL;
566- nsteps=-1;
567- analysis_type=-1;
568+ penalty = false;
569+ sid = -1;
570+ nodeid = -1;
571+ dof = -1;
572+ values = NULL;
573+ times = NULL;
574+ nsteps = -1;
575+ analysis_type = -1;
576 return;
577 }
578 /*}}}*/
579 SpcTransient::SpcTransient(int spc_sid,int spc_nodeid, int spc_dof,int spc_nsteps, IssmDouble* spc_times, IssmDouble* spc_values,int spc_analysis_type){/*{{{*/
580
581- sid=spc_sid;
582- nodeid=spc_nodeid;
583- dof=spc_dof;
584- nsteps=spc_nsteps;
585+ penalty= false;
586+ sid = spc_sid;
587+ nodeid = spc_nodeid;
588+ dof = spc_dof;
589+ nsteps = spc_nsteps;
590 if(spc_nsteps){
591- values=xNew<IssmDouble>(spc_nsteps);
592- times=xNew<IssmDouble>(spc_nsteps);
593+ values = xNew<IssmDouble>(spc_nsteps);
594+ times = xNew<IssmDouble>(spc_nsteps);
595 xMemCpy<IssmDouble>(values,spc_values,nsteps);
596 xMemCpy<IssmDouble>(times,spc_times,nsteps);
597 }
598@@ -68,7 +70,9 @@
599 this->Echo();
600 }
601 /*}}}*/
602-int SpcTransient::Id(void){ return sid; }/*{{{*/
603+int SpcTransient::Id(void){/*{{{*/
604+ return sid;
605+}
606 /*}}}*/
607 int SpcTransient::ObjectEnum(void){/*{{{*/
608
609@@ -100,7 +104,7 @@
610 /*Chase through nodes and find the node to which this SpcTransient applys: */
611 node=(Node*)nodes->GetObjectById(NULL,nodeid);
612
613- if(node){ //in case the spc is dealing with a node on another cpu
614+ if(!this->penalty && node){ //in case the spc is dealing with a node on another cpu
615
616 /*Retrieve time in parameters: */
617 parameters->FindParam(&time,TimeEnum);
618@@ -136,7 +140,63 @@
619 }
620 }
621 /*}}}*/
622+void SpcTransient::PenaltyDofAndValue(int* pdof,IssmDouble* pvalue,Nodes* nodes,Parameters* parameters){/*{{{*/
623
624+ if(!this->penalty) _error_("cannot return dof and value for non penalty constraint");
625+
626+ Node *node = NULL;
627+ IssmDouble time = 0.;
628+ int i,gdof;
629+ IssmDouble alpha = -1.;
630+ IssmDouble value;
631+ bool found = false;
632+
633+ /*Chase through nodes and find the node to which this SpcTransient applys: */
634+ node=(Node*)nodes->GetObjectById(NULL,nodeid);
635+
636+ if(node){ //in case the spc is dealing with a node on another cpu
637+
638+ /*Retrieve time in parameters: */
639+ parameters->FindParam(&time,TimeEnum);
640+
641+ /*Now, go fetch value for this time: */
642+ if (time<=times[0]){
643+ value=values[0];
644+ found=true;
645+ }
646+ else if (time>=times[nsteps-1]){
647+ value=values[nsteps-1];
648+ found=true;
649+ }
650+ else{
651+ for(i=0;i<nsteps-1;i++){
652+ if (times[i]<=time && time<times[i+1]){
653+ alpha=(time-times[i])/(times[i+1]-times[i]);
654+ value=(1-alpha)*values[i]+alpha*values[i+1];
655+ found=true;
656+ break;
657+ }
658+ }
659+ }
660+ if(!found)_error_("could not find time segment for constraint");
661+
662+ /*Get gdof */
663+ gdof = node->GetDof(dof,GsetEnum);
664+ if(xIsNan<IssmDouble>(value)){
665+ gdof = -1;
666+ }
667+ }
668+ else{
669+ value = NAN;
670+ gdof = -1;
671+ }
672+
673+ /*Assign output pointers*/
674+ *pdof = gdof;
675+ *pvalue = value;
676+}
677+/*}}}*/
678+
679 /*SpcTransient functions*/
680 int SpcTransient::GetDof(){/*{{{*/
681 return dof;
682Index: ../trunk-jpl/src/c
683===================================================================
684--- ../trunk-jpl/src/c (revision 18343)
685+++ ../trunk-jpl/src/c (revision 18344)
686
687Property changes on: ../trunk-jpl/src/c
688___________________________________________________________________
689Modified: svn:ignore
690## -4,6 +4,8 ##
691 *.o
692 *.obj
693 *.exe
694+issm
695+kriging
696 appscan.*
697 ouncemake*
698 list
Note: See TracBrowser for help on using the repository browser.