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

Contents of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 4 months ago) by jgs
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 40582 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26