source: issm/trunk-jpl/src/c/objects/Loads/Penpair.cpp@ 11332

Last change on this file since 11332 was 11332, checked in by Mathieu Morlighem, 13 years ago

Added Newton's method for MacAyeal
File size: 12.2 KB
Line 
1/*!\file Penpair.c
2 * \brief: implementation of the Penpair object
3 */
4
5/*Headers*/
6/*{{{1*/
7#ifdef HAVE_CONFIG_H
8#include <config.h>
9#else
10#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
11#endif
12
13#include <stdio.h>
14#include <string.h>
15#include "../objects.h"
16#include "../../EnumDefinitions/EnumDefinitions.h"
17#include "../../include/include.h"
18#include "../../shared/shared.h"
19/*}}}*/
20
21/*Element macros*/
22#define NUMVERTICES 2
23
24/*Penpair constructors and destructor*/
25/*FUNCTION Penpair::constructor {{{1*/
26Penpair::Penpair(){
27
28 this->hnodes=NULL;
29 this->nodes=NULL;
30 this->parameters=NULL;
31 return;
32}
33/*}}}1*/
34/*FUNCTION Penpair::creation {{{1*/
35Penpair::Penpair(int penpair_id, int* penpair_node_ids,int in_analysis_type){
36
37 this->id=penpair_id;
38 this->analysis_type=in_analysis_type;
39 this->hnodes=new Hook(penpair_node_ids,2);
40 this->parameters=NULL;
41 this->nodes=NULL;
42
43 return;
44}
45/*}}}1*/
46/*FUNCTION Penpair::destructor {{{1*/
47Penpair::~Penpair(){
48 delete hnodes;
49 return;
50}
51/*}}}1*/
52
53/*Object virtual functions definitions:*/
54/*FUNCTION Penpair::Echo {{{1*/
55void Penpair::Echo(void){
56
57 int i;
58
59 printf("Penpair:\n");
60 printf(" id: %i\n",id);
61 printf(" analysis_type: %s\n",EnumToStringx(analysis_type));
62 hnodes->Echo();
63
64 return;
65}
66/*}}}1*/
67/*FUNCTION Penpair::DeepEcho {{{1*/
68void Penpair::DeepEcho(void){
69
70 printf("Penpair:\n");
71 printf(" id: %i\n",id);
72 printf(" analysis_type: %s\n",EnumToStringx(analysis_type));
73 hnodes->DeepEcho();
74
75 return;
76}
77/*}}}1*/
78/*FUNCTION Penpair::Id {{{1*/
79int Penpair::Id(void){ return id; }
80/*}}}1*/
81/*FUNCTION Penpair::MyRank {{{1*/
82int Penpair::MyRank(void){
83 extern int my_rank;
84 return my_rank;
85}
86/*}}}1*/
87#ifdef _SERIAL_
88/*FUNCTION Penpair::Marshall {{{1*/
89void Penpair::Marshall(char** pmarshalled_dataset){
90
91 char* marshalled_dataset=NULL;
92 int enum_type=0;
93
94 /*recover marshalled_dataset: */
95 marshalled_dataset=*pmarshalled_dataset;
96
97 /*get enum type of Penpair: */
98 enum_type=PenpairEnum;
99
100 /*marshall enum: */
101 memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
102
103 /*marshall Penpair data: */
104 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
105 memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
106
107 /*Marshall hooks*/
108 hnodes->Marshall(&marshalled_dataset);
109
110 *pmarshalled_dataset=marshalled_dataset;
111 return;
112}
113/*}}}1*/
114/*FUNCTION Penpair::MarshallSize {{{1*/
115int Penpair::MarshallSize(){
116
117 return sizeof(id)+
118 +sizeof(analysis_type)
119 +hnodes->MarshallSize()
120 +sizeof(int); //sizeof(int) for enum type
121}
122/*}}}1*/
123/*FUNCTION Penpair::Demarshall {{{1*/
124void Penpair::Demarshall(char** pmarshalled_dataset){
125
126 int i;
127 char* marshalled_dataset=NULL;
128
129 /*recover marshalled_dataset: */
130 marshalled_dataset=*pmarshalled_dataset;
131
132 /*this time, no need to get enum type, the pointer directly points to the beginning of the
133 *object data (thanks to DataSet::Demarshall):*/
134 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
135 memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type);
136
137 /*demarshall hooks: */
138 hnodes=new Hook(); hnodes->Demarshall(&marshalled_dataset);
139
140 /*pointers are garbabe, until configuration is carried out: */
141 nodes=NULL;
142 parameters=NULL;
143
144 /*return: */
145 *pmarshalled_dataset=marshalled_dataset;
146 return;
147}
148/*}}}1*/
149#endif
150/*FUNCTION Penpair::ObjectEnum{{{1*/
151int Penpair::ObjectEnum(void){
152
153 return PenpairEnum;
154}
155/*}}}1*/
156/*FUNCTION Penpair::copy {{{1*/
157Object* Penpair::copy() {
158
159 Penpair* penpair=NULL;
160
161 penpair=new Penpair();
162
163 /*copy fields: */
164 penpair->id=this->id;
165 penpair->analysis_type=this->analysis_type;
166
167 /*now deal with hooks and objects: */
168 penpair->hnodes=(Hook*)this->hnodes->copy();
169 penpair->nodes =(Node**)penpair->hnodes->deliverp();
170
171 /*point parameters: */
172 penpair->parameters=this->parameters;
173
174 return penpair;
175
176}
177/*}}}*/
178
179/*Load virtual functions definitions:*/
180/*FUNCTION Penpair::Configure {{{1*/
181void Penpair::Configure(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
182
183 /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
184 * datasets, using internal ids and offsets hidden in hooks: */
185 hnodes->configure(nodesin);
186
187 /*Initialize hooked fields*/
188 this->nodes =(Node**)hnodes->deliverp();
189
190 /*point parameters to real dataset: */
191 this->parameters=parametersin;
192
193}
194/*}}}1*/
195/*FUNCTION Penpair::SetCurrentConfiguration {{{1*/
196void Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){
197
198}
199/*}}}1*/
200/*FUNCTION Penpair::CreateKMatrix {{{1*/
201void Penpair::CreateKMatrix(Mat Kff, Mat Kfs){
202 /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/
203 /*No loads applied, do nothing: */
204 return;
205
206}
207/*}}}1*/
208/*FUNCTION Penpair::CreatePVector {{{1*/
209void Penpair::CreatePVector(Vec pf){
210
211 /*No loads applied, do nothing: */
212 return;
213
214}
215/*}}}1*/
216/*FUNCTION Penpair::CreateJacobianMatrix{{{1*/
217void Penpair::CreateJacobianMatrix(Mat Jff){
218 this->CreateKMatrix(Jff,NULL);
219}
220/*}}}1*/
221/*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/
222void Penpair::PenaltyCreateKMatrix(Mat Kff, Mat Kfs,double kmax){
223
224 /*Retrieve parameters: */
225 ElementMatrix* Ke=NULL;
226 int analysis_type;
227 this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
228
229 switch(analysis_type){
230 case DiagnosticHorizAnalysisEnum:
231 Ke=PenaltyCreateKMatrixDiagnosticHoriz(kmax);
232 break;
233 case PrognosticAnalysisEnum:
234 Ke=PenaltyCreateKMatrixPrognostic(kmax);
235 break;
236 default:
237 _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
238 }
239
240 /*Add to global Vector*/
241 if(Ke){
242 Ke->AddToGlobal(Kff,Kfs);
243 delete Ke;
244 }
245}
246/*}}}1*/
247/*FUNCTION Penpair::PenaltyCreatePVector {{{1*/
248void Penpair::PenaltyCreatePVector(Vec pf,double kmax){
249 /*No loads applied, do nothing: */
250 return;
251}
252/*}}}1*/
253/*FUNCTION Penpair::PenaltyCreateJacobianMatrix{{{1*/
254void Penpair::PenaltyCreateJacobianMatrix(Mat Jff,double kmax){
255 this->PenaltyCreateKMatrix(Jff,NULL,kmax);
256}
257/*}}}1*/
258/*FUNCTION Penpair::InAnalysis{{{1*/
259bool Penpair::InAnalysis(int in_analysis_type){
260 if (in_analysis_type==this->analysis_type)return true;
261 else return false;
262}
263/*}}}*/
264
265/*Update virtual functions definitions:*/
266/*FUNCTION Penpair::InputUpdateFromConstant(double constant, int name) {{{1*/
267void Penpair::InputUpdateFromConstant(double constant, int name){
268 /*Nothing updated yet*/
269}
270/*}}}*/
271/*FUNCTION Penpair::InputUpdateFromConstant(int constant, int name) {{{1*/
272void Penpair::InputUpdateFromConstant(int constant, int name){
273 /*Nothing updated yet*/
274}
275/*}}}*/
276/*FUNCTION Penpair::InputUpdateFromConstant(bool constant, int name) {{{1*/
277void Penpair::InputUpdateFromConstant(bool constant, int name){
278 /*Nothing updated yet*/
279}
280/*}}}*/
281/*FUNCTION Penpair::InputUpdateFromVector(double* vector, int name, int type) {{{1*/
282void Penpair::InputUpdateFromVector(double* vector, int name, int type){
283 /*Nothing updated yet*/
284}
285/*}}}*/
286/*FUNCTION Penpair::InputUpdateFromVector(int* vector, int name, int type) {{{1*/
287void Penpair::InputUpdateFromVector(int* vector, int name, int type){
288 /*Nothing updated yet*/
289}
290/*}}}*/
291/*FUNCTION Penpair::InputUpdateFromVector(bool* vector, int name, int type) {{{1*/
292void Penpair::InputUpdateFromVector(bool* vector, int name, int type){
293 /*Nothing updated yet*/
294}
295/*}}}*/
296
297/*Penpair management:*/
298/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticHoriz{{{1*/
299ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticHoriz(double kmax){
300
301 int approximation0=nodes[0]->GetApproximation();
302 int approximation1=nodes[1]->GetApproximation();
303
304 switch(approximation0){
305 case MacAyealApproximationEnum:
306 switch(approximation1){
307 case MacAyealApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax);
308 case PattynApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax);
309 default: _error_("not supported yet");
310 }
311 case PattynApproximationEnum:
312 switch(approximation1){
313 case MacAyealApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax);
314 case PattynApproximationEnum: return PenaltyCreateKMatrixDiagnosticMacAyealPattyn(kmax);
315 default: _error_("not supported yet");
316 }
317 case StokesApproximationEnum:
318 switch(approximation1){
319 case StokesApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax);
320 case NoneApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax);
321 default: _error_("not supported yet");
322 }
323 case NoneApproximationEnum:
324 switch(approximation1){
325 case StokesApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax);
326 case NoneApproximationEnum: return PenaltyCreateKMatrixDiagnosticStokes(kmax);
327 default: _error_("not supported yet");
328 }
329 default: _error_("not supported yet");
330 }
331}
332/*}}}1*/
333/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn {{{1*/
334ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticMacAyealPattyn(double kmax){
335
336 const int numdof=NUMVERTICES*NDOF2;
337 double penalty_offset;
338
339 /*Initialize Element vector and return if necessary*/
340 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
341
342 /*recover parameters: */
343 parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);
344
345 //Create elementary matrix: add penalty to
346 Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_offset);
347 Ke->values[0*numdof+2]=-kmax*pow((double)10.0,penalty_offset);
348 Ke->values[2*numdof+0]=-kmax*pow((double)10.0,penalty_offset);
349 Ke->values[2*numdof+2]=+kmax*pow((double)10.0,penalty_offset);
350
351 Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_offset);
352 Ke->values[1*numdof+3]=-kmax*pow((double)10.0,penalty_offset);
353 Ke->values[3*numdof+1]=-kmax*pow((double)10.0,penalty_offset);
354 Ke->values[3*numdof+3]=+kmax*pow((double)10.0,penalty_offset);
355
356 /*Clean up and return*/
357 return Ke;
358}
359/*}}}1*/
360/*FUNCTION Penpair::PenaltyCreateKMatrixDiagnosticStokes {{{1*/
361ElementMatrix* Penpair::PenaltyCreateKMatrixDiagnosticStokes(double kmax){
362
363 const int numdof=NUMVERTICES*NDOF4;
364 double penalty_offset;
365
366 /*Initialize Element vector and return if necessary*/
367 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
368
369 /*recover parameters: */
370 parameters->FindParam(&penalty_offset,DiagnosticPenaltyFactorEnum);
371
372 //Create elementary matrix: add penalty to
373 Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_offset);
374 Ke->values[0*numdof+4]=-kmax*pow((double)10.0,penalty_offset);
375 Ke->values[4*numdof+0]=-kmax*pow((double)10.0,penalty_offset);
376 Ke->values[4*numdof+4]=+kmax*pow((double)10.0,penalty_offset);
377
378 Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_offset);
379 Ke->values[1*numdof+5]=-kmax*pow((double)10.0,penalty_offset);
380 Ke->values[5*numdof+1]=-kmax*pow((double)10.0,penalty_offset);
381 Ke->values[5*numdof+5]=+kmax*pow((double)10.0,penalty_offset);
382
383 Ke->values[2*numdof+2]=+kmax*pow((double)10.0,penalty_offset);
384 Ke->values[2*numdof+6]=-kmax*pow((double)10.0,penalty_offset);
385 Ke->values[6*numdof+2]=-kmax*pow((double)10.0,penalty_offset);
386 Ke->values[6*numdof+6]=+kmax*pow((double)10.0,penalty_offset);
387
388 Ke->values[3*numdof+3]=+kmax*pow((double)10.0,penalty_offset);
389 Ke->values[3*numdof+7]=-kmax*pow((double)10.0,penalty_offset);
390 Ke->values[7*numdof+3]=-kmax*pow((double)10.0,penalty_offset);
391 Ke->values[7*numdof+7]=+kmax*pow((double)10.0,penalty_offset);
392
393 /*Clean up and return*/
394 return Ke;
395}
396/*}}}1*/
397/*FUNCTION Penpair::PenaltyCreateKMatrixPrognostic {{{1*/
398ElementMatrix* Penpair::PenaltyCreateKMatrixPrognostic(double kmax){
399
400 const int numdof=NUMVERTICES*NDOF1;
401 double penalty_factor;
402
403 /*Initialize Element vector and return if necessary*/
404 ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters);
405
406 /*recover parameters: */
407 parameters->FindParam(&penalty_factor,PrognosticPenaltyFactorEnum);
408
409 //Create elementary matrix: add penalty to
410 Ke->values[0*numdof+0]=+kmax*pow((double)10.0,penalty_factor);
411 Ke->values[0*numdof+1]=-kmax*pow((double)10.0,penalty_factor);
412 Ke->values[1*numdof+0]=-kmax*pow((double)10.0,penalty_factor);
413 Ke->values[1*numdof+1]=+kmax*pow((double)10.0,penalty_factor);
414
415 /*Clean up and return*/
416 return Ke;
417}
418/*}}}1*/
Note: See TracBrowser for help on using the repository browser.