/[escript]/trunk/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 472 - (show annotations)
Fri Jan 27 01:50:59 2006 UTC (12 years, 11 months ago) by jgs
File size: 39959 byte(s)
rationalise all #includes

1 // $Id$
2 /*
3 ******************************************************************************
4 * *
5 * COPYRIGHT ACcESS 2004 - All Rights Reserved *
6 * *
7 * This software is the property of ACcESS. No part of this code *
8 * may be copied in any form or by any means without the expressed written *
9 * consent of ACcESS. Copying, use or modification of this software *
10 * by any unauthorised person is illegal unless that person has a software *
11 * license agreement with ACcESS. *
12 * *
13 ******************************************************************************
14 */
15
16 #include "MeshAdapter.h"
17
18 using namespace std;
19 using namespace escript;
20
21 namespace finley {
22
23 //
24 // define the static constants
25 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
26 const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
27 const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
28 const int MeshAdapter::Nodes=FINLEY_NODES;
29 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
30 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
31 const int MeshAdapter::Points=FINLEY_POINTS;
32 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
33 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
34
35 MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
36 {
37 setFunctionSpaceTypeNames();
38 //
39 // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer
40 // for us.
41 m_finleyMesh.reset(finleyMesh,null_deleter());
42 }
43
44 //
45 // The copy constructor should just increment the use count
46 MeshAdapter::MeshAdapter(const MeshAdapter& in):
47 m_finleyMesh(in.m_finleyMesh)
48 {
49 setFunctionSpaceTypeNames();
50 }
51
52 MeshAdapter::~MeshAdapter()
53 {
54 //
55 // I hope the case for the pointer being zero has been taken care of.
56 // cout << "In MeshAdapter destructor." << endl;
57 if (m_finleyMesh.unique()) {
58 // cout << "Calling dealloc." << endl;
59 Finley_Mesh_dealloc(m_finleyMesh.get());
60 // cout << "Finished dealloc." << endl;
61 }
62 }
63
64 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
65 return m_finleyMesh.get();
66 }
67
68 void MeshAdapter::write(const std::string& fileName) const
69 {
70 char fName[fileName.size()+1];
71 strcpy(fName,fileName.c_str());
72 Finley_Mesh_write(m_finleyMesh.get(),fName);
73 checkFinleyError();
74 }
75
76 // void MeshAdapter::getTagList(int functionSpaceType,
77 // int* numTags) const
78 // {
79 // Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);
80 // return;
81 // }
82
83 string MeshAdapter::getDescription() const
84 {
85 return "FinleyMesh";
86 }
87
88 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
89 {
90 FunctionSpaceNamesMapType::iterator loc;
91 loc=m_functionSpaceTypeNames.find(functionSpaceType);
92 if (loc==m_functionSpaceTypeNames.end()) {
93 return "Invalid function space type code.";
94 } else {
95 return loc->second;
96 }
97 }
98
99 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
100 {
101 FunctionSpaceNamesMapType::iterator loc;
102 loc=m_functionSpaceTypeNames.find(functionSpaceType);
103 return (loc!=m_functionSpaceTypeNames.end());
104 }
105
106 void MeshAdapter::setFunctionSpaceTypeNames()
107 {
108 m_functionSpaceTypeNames.insert
109 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
110 m_functionSpaceTypeNames.insert
111 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
112 m_functionSpaceTypeNames.insert
113 (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
114 m_functionSpaceTypeNames.insert
115 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
116 m_functionSpaceTypeNames.insert
117 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
118 m_functionSpaceTypeNames.insert
119 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
120 m_functionSpaceTypeNames.insert
121 (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
122 m_functionSpaceTypeNames.insert
123 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
124 }
125
126 int MeshAdapter::getContinuousFunctionCode() const
127 {
128 return Nodes;
129 }
130
131 int MeshAdapter::getFunctionCode() const
132 {
133 return Elements;
134 }
135
136 int MeshAdapter::getFunctionOnBoundaryCode() const
137 {
138 return FaceElements;
139 }
140
141 int MeshAdapter::getFunctionOnContactZeroCode() const
142 {
143 return ContactElementsZero;
144 }
145
146 int MeshAdapter::getFunctionOnContactOneCode() const
147 {
148 return ContactElementsOne;
149 }
150
151 int MeshAdapter::getSolutionCode() const
152 {
153 return DegreesOfFreedom;
154 }
155
156 int MeshAdapter::getReducedSolutionCode() const
157 {
158 return ReducedDegreesOfFreedom;
159 }
160
161 int MeshAdapter::getDiracDeltaFunctionCode() const
162 {
163 return Points;
164 }
165
166 //
167 // returns a pointer to the tag list of samples of functionSpaceType
168 //
169 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,
170 int* numTags) const
171 {
172 *tagList=NULL;
173 *numTags=0;
174 Finley_Mesh* mesh=m_finleyMesh.get();
175 switch (functionSpaceType) {
176 case(Nodes):
177 if (mesh->Nodes!=NULL) {
178 *tagList=mesh->Nodes->Tag;
179 *numTags=mesh->Nodes->numNodes;
180 }
181 break;
182 case(Elements):
183 if (mesh->Elements!=NULL) {
184 *tagList=mesh->Elements->Tag;
185 *numTags=mesh->Elements->numElements;
186 }
187 break;
188 case(FaceElements):
189 if (mesh->FaceElements!=NULL) {
190 *tagList=mesh->FaceElements->Tag;
191 *numTags=mesh->FaceElements->numElements;
192 }
193 break;
194 case(Points):
195 if (mesh->Points!=NULL) {
196 *tagList=mesh->Points->Tag;
197 *numTags=mesh->Points->numElements;
198 }
199 break;
200 case(ContactElementsZero):
201 if (mesh->ContactElements!=NULL) {
202 *tagList=mesh->ContactElements->Tag;
203 *numTags=mesh->ContactElements->numElements;
204 }
205 break;
206 case(ContactElementsOne):
207 if (mesh->ContactElements!=NULL) {
208 *tagList=mesh->ContactElements->Tag;
209 *numTags=mesh->ContactElements->numElements;
210 }
211 break;
212 case(DegreesOfFreedom):
213 if (mesh->Nodes!=NULL) {
214 *tagList=NULL;
215 *numTags=0;
216 }
217 break;
218 case(ReducedDegreesOfFreedom):
219 if (mesh->Nodes!=NULL) {
220 *tagList=NULL;
221 *numTags=0;
222 }
223 break;
224 default:
225 stringstream temp;
226 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
227 throw FinleyAdapterException(temp.str());
228 break;
229 }
230 if (*tagList==NULL) {
231 stringstream temp;
232 temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();
233 throw FinleyAdapterException(temp.str());
234 }
235 return;
236 }
237
238 //
239 // returns a pointer to the reference no list of samples of functionSpaceType
240 //
241 void MeshAdapter::getReferenceNoList(int functionSpaceType, int** referenceNoList,
242 int* numReferenceNo) const
243 {
244 *referenceNoList=NULL;
245 *numReferenceNo=0;
246 Finley_Mesh* mesh=m_finleyMesh.get();
247 switch (functionSpaceType) {
248 case(Nodes):
249 if (mesh->Nodes!=NULL) {
250 *referenceNoList=mesh->Nodes->Id;
251 *numReferenceNo=mesh->Nodes->numNodes;
252 }
253 break;
254 case(Elements):
255 if (mesh->Elements!=NULL) {
256 *referenceNoList=mesh->Elements->Id;
257 *numReferenceNo=mesh->Elements->numElements;
258 }
259 break;
260 case(FaceElements):
261 if (mesh->FaceElements!=NULL) {
262 *referenceNoList=mesh->FaceElements->Id;
263 *numReferenceNo=mesh->FaceElements->numElements;
264 }
265 break;
266 case(Points):
267 if (mesh->Points!=NULL) {
268 *referenceNoList=mesh->Points->Id;
269 *numReferenceNo=mesh->Points->numElements;
270 }
271 break;
272 case(ContactElementsZero):
273 if (mesh->ContactElements!=NULL) {
274 *referenceNoList=mesh->ContactElements->Id;
275 *numReferenceNo=mesh->ContactElements->numElements;
276 }
277 break;
278 case(ContactElementsOne):
279 if (mesh->ContactElements!=NULL) {
280 *referenceNoList=mesh->ContactElements->Id;
281 *numReferenceNo=mesh->ContactElements->numElements;
282 }
283 break;
284 case(DegreesOfFreedom):
285 if (mesh->Nodes!=NULL) {
286 *referenceNoList=NULL;
287 *numReferenceNo=0;
288 }
289 break;
290 case(ReducedDegreesOfFreedom):
291 if (mesh->Nodes!=NULL) {
292 *referenceNoList=NULL;
293 *numReferenceNo=0;
294 }
295 break;
296 default:
297 stringstream temp;
298 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
299 throw FinleyAdapterException(temp.str());
300 break;
301 }
302 if (*referenceNoList==NULL) {
303 stringstream temp;
304 temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();
305 throw FinleyAdapterException(temp.str());
306 }
307 return;
308 }
309
310 //
311 // return the spatial dimension of the Mesh:
312 //
313 int MeshAdapter::getDim() const
314 {
315 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
316 checkFinleyError();
317 return numDim;
318 }
319
320 //
321 // return the number of data points per sample and the number of samples
322 // needed to represent data on a parts of the mesh.
323 //
324 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
325 {
326 int numDataPointsPerSample=0;
327 int numSamples=0;
328 Finley_Mesh* mesh=m_finleyMesh.get();
329 switch (functionSpaceCode) {
330 case(Nodes):
331 numDataPointsPerSample=1;
332 if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
333 break;
334 case(Elements):
335 if (mesh->Elements!=NULL) {
336 numSamples=mesh->Elements->numElements;
337 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
338 }
339 break;
340 case(FaceElements):
341 if (mesh->FaceElements!=NULL) {
342 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
343 numSamples=mesh->FaceElements->numElements;
344 }
345 break;
346 case(Points):
347 if (mesh->Points!=NULL) {
348 numDataPointsPerSample=1;
349 numSamples=mesh->Points->numElements;
350 }
351 break;
352 case(ContactElementsZero):
353 if (mesh->ContactElements!=NULL) {
354 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
355 numSamples=mesh->ContactElements->numElements;
356 }
357 break;
358 case(ContactElementsOne):
359 if (mesh->ContactElements!=NULL) {
360 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
361 numSamples=mesh->ContactElements->numElements;
362 }
363 break;
364 case(DegreesOfFreedom):
365 if (mesh->Nodes!=NULL) {
366 numDataPointsPerSample=1;
367 numSamples=mesh->Nodes->numDegreesOfFreedom;
368 }
369 break;
370 case(ReducedDegreesOfFreedom):
371 if (mesh->Nodes!=NULL) {
372 numDataPointsPerSample=1;
373 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
374 }
375 break;
376 default:
377 stringstream temp;
378 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
379 throw FinleyAdapterException(temp.str());
380 break;
381 }
382 return pair<int,int>(numDataPointsPerSample,numSamples);
383 }
384
385 //
386 // adds linear PDE of second order into a given stiffness matrix and right hand side:
387 //
388 void MeshAdapter::addPDEToSystem(
389 SystemMatrixAdapter& mat, Data& rhs,
390 const Data& A, const Data& B, const Data& C,const Data& D,const Data& X,const Data& Y,
391 const Data& d, const Data& y,
392 const Data& d_contact,const Data& y_contact) const
393 {
394 Finley_Mesh* mesh=m_finleyMesh.get();
395 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),
396 &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
397 checkFinleyError();
398 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
399 mat.getPaso_SystemMatrix(),
400 &(rhs.getDataC()),
401 &(d.getDataC()),&(y.getDataC()),
402 Finley_Assemble_handelShapeMissMatch_Mean_out);
403 checkFinleyError();
404 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,
405 mat.getPaso_SystemMatrix(),
406 &(rhs.getDataC()),
407 &(d_contact.getDataC()),
408 &(y_contact.getDataC()),
409 Finley_Assemble_handelShapeMissMatch_Step_out);
410 checkFinleyError();
411 }
412
413 //
414 // adds linear PDE of second order into the right hand side only
415 //
416 void MeshAdapter::addPDEToRHS( Data& rhs,
417 const Data& X,const Data& Y, const Data& y, const Data& y_contact) const
418 {
419 Finley_Mesh* mesh=m_finleyMesh.get();
420
421 // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));
422 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));
423 checkFinleyError();
424
425 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
426 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
427
428 checkFinleyError();
429 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
430 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
431 checkFinleyError();
432 }
433
434 //
435 // interpolates data between different function spaces:
436 //
437 void MeshAdapter::interpolateOnDomain(Data& target,const Data& in) const
438 {
439 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
440 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
441 if (inDomain!=*this)
442 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
443 if (targetDomain!=*this)
444 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
445
446 Finley_Mesh* mesh=m_finleyMesh.get();
447 switch(in.getFunctionSpace().getTypeCode()) {
448 case(Nodes):
449 switch(target.getFunctionSpace().getTypeCode()) {
450 case(Nodes):
451 case(ReducedDegreesOfFreedom):
452 case(DegreesOfFreedom):
453 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
454 break;
455 case(Elements):
456 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
457 break;
458 case(FaceElements):
459 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
460 break;
461 case(Points):
462 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
463 break;
464 case(ContactElementsZero):
465 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
466 break;
467 case(ContactElementsOne):
468 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
469 break;
470 default:
471 stringstream temp;
472 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
473 throw FinleyAdapterException(temp.str());
474 break;
475 }
476 break;
477 case(Elements):
478 if (target.getFunctionSpace().getTypeCode()==Elements) {
479 Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
480 } else {
481 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
482 }
483 break;
484 case(FaceElements):
485 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
486 Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
487 } else {
488 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
489 break;
490 }
491 case(Points):
492 if (target.getFunctionSpace().getTypeCode()==Points) {
493 Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
494 } else {
495 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
496 break;
497 }
498 break;
499 case(ContactElementsZero):
500 case(ContactElementsOne):
501 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
502 Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
503 } else {
504 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
505 break;
506 }
507 break;
508 case(DegreesOfFreedom):
509 switch(target.getFunctionSpace().getTypeCode()) {
510 case(ReducedDegreesOfFreedom):
511 case(DegreesOfFreedom):
512 case(Nodes):
513 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
514 break;
515 case(Elements):
516 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
517 break;
518 case(FaceElements):
519 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
520 break;
521 case(Points):
522 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
523 break;
524 case(ContactElementsZero):
525 case(ContactElementsOne):
526 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
527 break;
528 default:
529 stringstream temp;
530 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
531 throw FinleyAdapterException(temp.str());
532 break;
533 }
534 break;
535 case(ReducedDegreesOfFreedom):
536 switch(target.getFunctionSpace().getTypeCode()) {
537 case(ReducedDegreesOfFreedom):
538 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
539 break;
540 case(Elements):
541 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
542 break;
543 case(FaceElements):
544 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
545 break;
546 case(Points):
547 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
548 break;
549 case(ContactElementsZero):
550 case(ContactElementsOne):
551 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
552 break;
553 case(Nodes):
554 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
555 break;
556 case(DegreesOfFreedom):
557 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
558 break;
559 default:
560 stringstream temp;
561 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
562 throw FinleyAdapterException(temp.str());
563 break;
564 }
565 break;
566 default:
567 stringstream temp;
568 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
569 throw FinleyAdapterException(temp.str());
570 break;
571 }
572 checkFinleyError();
573 }
574
575 //
576 // copies the locations of sample points into x:
577 //
578 void MeshAdapter::setToX(Data& arg) const
579 {
580 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
581 if (argDomain!=*this)
582 throw FinleyAdapterException("Error - Illegal domain of data point locations");
583 Finley_Mesh* mesh=m_finleyMesh.get();
584 // in case of values node coordinates we can do the job directly:
585 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
586 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
587 } else {
588 Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
589 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
590 // this is then interpolated onto arg:
591 interpolateOnDomain(arg,tmp_data);
592 }
593 checkFinleyError();
594 }
595
596 //
597 // return the normal vectors at the location of data points as a Data object:
598 //
599 void MeshAdapter::setToNormal(Data& normal) const
600 {
601 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
602 if (normalDomain!=*this)
603 throw FinleyAdapterException("Error - Illegal domain of normal locations");
604 Finley_Mesh* mesh=m_finleyMesh.get();
605 switch(normal.getFunctionSpace().getTypeCode()) {
606 case(Nodes):
607 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
608 break;
609 case(Elements):
610 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
611 break;
612 case (FaceElements):
613 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
614 break;
615 case(Points):
616 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
617 break;
618 case (ContactElementsOne):
619 case (ContactElementsZero):
620 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
621 break;
622 case(DegreesOfFreedom):
623 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
624 break;
625 case(ReducedDegreesOfFreedom):
626 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
627 break;
628 default:
629 stringstream temp;
630 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
631 throw FinleyAdapterException(temp.str());
632 break;
633 }
634 checkFinleyError();
635 }
636
637 //
638 // interpolates data to other domain:
639 //
640 void MeshAdapter::interpolateACross(Data& target,const Data& source) const
641 {
642 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
643 if (targetDomain!=*this)
644 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
645
646 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
647 }
648
649 //
650 // calculates the integral of a function defined of arg:
651 //
652 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const Data& arg) const
653 {
654 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
655 if (argDomain!=*this)
656 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
657
658 Finley_Mesh* mesh=m_finleyMesh.get();
659 switch(arg.getFunctionSpace().getTypeCode()) {
660 case(Nodes):
661 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
662 break;
663 case(Elements):
664 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
665 break;
666 case(FaceElements):
667 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
668 break;
669 case(Points):
670 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
671 break;
672 case(ContactElementsZero):
673 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
674 break;
675 case(ContactElementsOne):
676 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
677 break;
678 case(DegreesOfFreedom):
679 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
680 break;
681 case(ReducedDegreesOfFreedom):
682 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
683 break;
684 default:
685 stringstream temp;
686 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
687 throw FinleyAdapterException(temp.str());
688 break;
689 }
690 checkFinleyError();
691 }
692
693 //
694 // calculates the gradient of arg:
695 //
696 void MeshAdapter::setToGradient(Data& grad,const Data& arg) const
697 {
698 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
699 if (argDomain!=*this)
700 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
701 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
702 if (gradDomain!=*this)
703 throw FinleyAdapterException("Error - Illegal domain of gradient");
704
705 Finley_Mesh* mesh=m_finleyMesh.get();
706 switch(grad.getFunctionSpace().getTypeCode()) {
707 case(Nodes):
708 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
709 break;
710 case(Elements):
711 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));
712 break;
713 case(FaceElements):
714 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));
715 break;
716 case(Points):
717 throw FinleyAdapterException("Error - Gradient at points is not supported.");
718 break;
719 case(ContactElementsZero):
720 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
721 break;
722 case(ContactElementsOne):
723 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
724 break;
725 case(DegreesOfFreedom):
726 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
727 break;
728 case(ReducedDegreesOfFreedom):
729 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
730 break;
731 default:
732 stringstream temp;
733 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
734 throw FinleyAdapterException(temp.str());
735 break;
736 }
737 checkFinleyError();
738 }
739
740 //
741 // returns the size of elements:
742 //
743 void MeshAdapter::setToSize(Data& size) const
744 {
745 Finley_Mesh* mesh=m_finleyMesh.get();
746 escriptDataC tmp=size.getDataC();
747 switch(size.getFunctionSpace().getTypeCode()) {
748 case(Nodes):
749 throw FinleyAdapterException("Error - Size of nodes is not supported.");
750 break;
751 case(Elements):
752 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
753 break;
754 case(FaceElements):
755 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
756 break;
757 case(Points):
758 throw FinleyAdapterException("Error - Size of point elements is not supported.");
759 break;
760 case(ContactElementsZero):
761 case(ContactElementsOne):
762 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
763 break;
764 case(DegreesOfFreedom):
765 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
766 break;
767 case(ReducedDegreesOfFreedom):
768 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
769 break;
770 default:
771 stringstream temp;
772 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
773 throw FinleyAdapterException(temp.str());
774 break;
775 }
776 checkFinleyError();
777 }
778
779 // sets the location of nodes:
780 void MeshAdapter::setNewX(const Data& new_x)
781 {
782 Finley_Mesh* mesh=m_finleyMesh.get();
783 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
784 if (newDomain!=*this)
785 throw FinleyAdapterException("Error - Illegal domain of new point locations");
786 Finley_NodeFile_setCoordinates(mesh->Nodes,&(new_x.getDataC()));
787 checkFinleyError();
788 }
789
790 // saves a data array in openDX format:
791 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
792 {
793 int MAX_namelength=256;
794 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
795 char names[num_data][MAX_namelength];
796 char* c_names[num_data];
797 escriptDataC data[num_data];
798 escriptDataC* ptr_data[num_data];
799
800 boost::python::list keys=arg.keys();
801 for (int i=0;i<num_data;++i) {
802 Data& d=boost::python::extract<Data&>(arg[keys[i]]);
803 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
804 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
805 data[i]=d.getDataC();
806 ptr_data[i]=&(data[i]);
807 std::string n=boost::python::extract<std::string>(keys[i]);
808 c_names[i]=&(names[i][0]);
809 if (n.length()>MAX_namelength-1) {
810 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
811 c_names[i][MAX_namelength-1]='\0';
812 } else {
813 strcpy(c_names[i],n.c_str());
814 }
815 }
816 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
817 checkFinleyError();
818 return;
819 }
820
821 // saves a data array in openVTK format:
822 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
823 {
824 int MAX_namelength=256;
825 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
826 char names[num_data][MAX_namelength];
827 char* c_names[num_data];
828 escriptDataC data[num_data];
829 escriptDataC* ptr_data[num_data];
830
831 boost::python::list keys=arg.keys();
832 for (int i=0;i<num_data;++i) {
833 Data& d=boost::python::extract<Data&>(arg[keys[i]]);
834 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
835 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
836 data[i]=d.getDataC();
837 ptr_data[i]=&(data[i]);
838 std::string n=boost::python::extract<std::string>(keys[i]);
839 c_names[i]=&(names[i][0]);
840 if (n.length()>MAX_namelength-1) {
841 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
842 c_names[i][MAX_namelength-1]='\0';
843 } else {
844 strcpy(c_names[i],n.c_str());
845 }
846 }
847 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
848 checkFinleyError();
849 return;
850 }
851
852
853 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
854 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
855 const int row_blocksize,
856 const escript::FunctionSpace& row_functionspace,
857 const int column_blocksize,
858 const escript::FunctionSpace& column_functionspace,
859 const int type) const
860 {
861 int reduceRowOrder=0;
862 int reduceColOrder=0;
863 // is the domain right?
864 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
865 if (row_domain!=*this)
866 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
867 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
868 if (col_domain!=*this)
869 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
870 // is the function space type right
871 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
872 reduceRowOrder=0;
873 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
874 reduceRowOrder=1;
875 } else {
876 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
877 }
878 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
879 reduceColOrder=0;
880 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
881 reduceColOrder=1;
882 } else {
883 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
884 }
885 // generate matrix:
886
887 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
888 checkFinleyError();
889 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
890 checkPasoError();
891 Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
892 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
893 }
894
895 //
896 // vtkObject MeshAdapter::createVtkObject() const
897 // TODO:
898 //
899
900 //
901 // returns true if data at the atom_type is considered as being cell centered:
902 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
903 {
904 switch(functionSpaceCode) {
905 case(Nodes):
906 case(DegreesOfFreedom):
907 case(ReducedDegreesOfFreedom):
908 return false;
909 break;
910 case(Elements):
911 case(FaceElements):
912 case(Points):
913 case(ContactElementsZero):
914 case(ContactElementsOne):
915 return true;
916 break;
917 default:
918 stringstream temp;
919 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
920 throw FinleyAdapterException(temp.str());
921 break;
922 }
923 checkFinleyError();
924 return false;
925 }
926
927 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
928 {
929 switch(functionSpaceType_source) {
930 case(Nodes):
931 switch(functionSpaceType_target) {
932 case(Nodes):
933 case(ReducedDegreesOfFreedom):
934 case(DegreesOfFreedom):
935 case(Elements):
936 case(FaceElements):
937 case(Points):
938 case(ContactElementsZero):
939 case(ContactElementsOne):
940 return true;
941 default:
942 stringstream temp;
943 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
944 throw FinleyAdapterException(temp.str());
945 }
946 break;
947 case(Elements):
948 if (functionSpaceType_target==Elements) {
949 return true;
950 } else {
951 return false;
952 }
953 case(FaceElements):
954 if (functionSpaceType_target==FaceElements) {
955 return true;
956 } else {
957 return false;
958 }
959 case(Points):
960 if (functionSpaceType_target==Points) {
961 return true;
962 } else {
963 return false;
964 }
965 case(ContactElementsZero):
966 case(ContactElementsOne):
967 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
968 return true;
969 } else {
970 return false;
971 }
972 case(DegreesOfFreedom):
973 switch(functionSpaceType_target) {
974 case(ReducedDegreesOfFreedom):
975 case(DegreesOfFreedom):
976 case(Nodes):
977 case(Elements):
978 case(FaceElements):
979 case(Points):
980 case(ContactElementsZero):
981 case(ContactElementsOne):
982 return true;
983 default:
984 stringstream temp;
985 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
986 throw FinleyAdapterException(temp.str());
987 }
988 break;
989 case(ReducedDegreesOfFreedom):
990 switch(functionSpaceType_target) {
991 case(ReducedDegreesOfFreedom):
992 case(Elements):
993 case(FaceElements):
994 case(Points):
995 case(ContactElementsZero):
996 case(ContactElementsOne):
997 return true;
998 case(Nodes):
999 case(DegreesOfFreedom):
1000 return false;
1001 default:
1002 stringstream temp;
1003 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1004 throw FinleyAdapterException(temp.str());
1005 }
1006 break;
1007 default:
1008 stringstream temp;
1009 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1010 throw FinleyAdapterException(temp.str());
1011 break;
1012 }
1013 checkFinleyError();
1014 return false;
1015 }
1016
1017 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1018 {
1019 return false;
1020 }
1021
1022 bool MeshAdapter::operator==(const AbstractDomain& other) const
1023 {
1024 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1025 if (temp!=0) {
1026 return (m_finleyMesh==temp->m_finleyMesh);
1027 } else {
1028 return false;
1029 }
1030 }
1031
1032 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1033 {
1034 return !(operator==(other));
1035 }
1036
1037 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1038 {
1039 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1040 checkPasoError();
1041 return out;
1042 }
1043
1044 Data MeshAdapter::getX() const
1045 {
1046 return continuousFunction(asAbstractContinuousDomain()).getX();
1047 }
1048
1049 Data MeshAdapter::getNormal() const
1050 {
1051 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1052 }
1053
1054 Data MeshAdapter::getSize() const
1055 {
1056 return function(asAbstractContinuousDomain()).getSize();
1057 }
1058
1059 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1060 {
1061 int* tagList;
1062 int numTags;
1063 getTagList(functionSpaceType, &tagList, &numTags);
1064 return tagList[sampleNo];
1065 }
1066
1067 int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1068 {
1069 int* referenceNoList;
1070 int numReferenceNo;
1071 getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);
1072 return referenceNoList[sampleNo];
1073 }
1074
1075 } // end of namespace

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26