/[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 1137 - (show annotations)
Thu May 10 08:11:31 2007 UTC (12 years, 1 month ago) by gross
File size: 56150 byte(s)
This version passes the tests on windows except for 

   * vtk
   * netCDF

The version needs to be tested on altix and linux
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 #ifdef PASO_MPI
17 #include <mpi.h>
18 #endif
19 #include "MeshAdapter.h"
20
21 #include "escript/Data.h"
22 #include "escript/DataFactory.h"
23
24 using namespace std;
25 using namespace escript;
26
27 namespace finley {
28
29 //
30 // define the static constants
31 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
32 const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
33 const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
34 const int MeshAdapter::Nodes=FINLEY_NODES;
35 const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
36 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
37 const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
38 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
39 const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
40 const int MeshAdapter::Points=FINLEY_POINTS;
41 const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
42 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
43 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
44 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
45
46 MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
47 {
48 setFunctionSpaceTypeNames();
49 //
50 // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer
51 // for us.
52 m_finleyMesh.reset(finleyMesh,null_deleter());
53 }
54
55 //
56 // The copy constructor should just increment the use count
57 MeshAdapter::MeshAdapter(const MeshAdapter& in):
58 m_finleyMesh(in.m_finleyMesh)
59 {
60 setFunctionSpaceTypeNames();
61 }
62
63 MeshAdapter::~MeshAdapter()
64 {
65 //
66 // I hope the case for the pointer being zero has been taken care of.
67 // cout << "In MeshAdapter destructor." << endl;
68 if (m_finleyMesh.unique()) {
69 // cout << "Calling dealloc." << endl;
70 Finley_Mesh_dealloc(m_finleyMesh.get());
71 // cout << "Finished dealloc." << endl;
72 }
73 }
74
75 Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
76 return m_finleyMesh.get();
77 }
78
79 void MeshAdapter::write(const std::string& fileName) const
80 {
81 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
82 strcpy(fName,fileName.c_str());
83 Finley_Mesh_write(m_finleyMesh.get(),fName);
84 checkFinleyError();
85 TMPMEMFREE(fName);
86 }
87
88 string MeshAdapter::getDescription() const
89 {
90 return string("FinleyMesh");
91 }
92
93 string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
94 {
95 FunctionSpaceNamesMapType::iterator loc;
96 loc=m_functionSpaceTypeNames.find(functionSpaceType);
97 if (loc==m_functionSpaceTypeNames.end()) {
98 return "Invalid function space type code.";
99 } else {
100 return string(loc->second);
101 }
102 }
103
104 bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
105 {
106 FunctionSpaceNamesMapType::iterator loc;
107 loc=m_functionSpaceTypeNames.find(functionSpaceType);
108 return (loc!=m_functionSpaceTypeNames.end());
109 }
110
111 void MeshAdapter::setFunctionSpaceTypeNames()
112 {
113 m_functionSpaceTypeNames.insert
114 (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
115 m_functionSpaceTypeNames.insert
116 (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
117 m_functionSpaceTypeNames.insert
118 (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
119 m_functionSpaceTypeNames.insert
120 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
121 m_functionSpaceTypeNames.insert
122 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
123 m_functionSpaceTypeNames.insert
124 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
125 m_functionSpaceTypeNames.insert
126 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
127 m_functionSpaceTypeNames.insert
128 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
129 m_functionSpaceTypeNames.insert
130 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
131 m_functionSpaceTypeNames.insert
132 (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
133 m_functionSpaceTypeNames.insert
134 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
135 m_functionSpaceTypeNames.insert
136 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
137 m_functionSpaceTypeNames.insert
138 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
139 }
140
141 int MeshAdapter::getContinuousFunctionCode() const
142 {
143 return Nodes;
144 }
145 int MeshAdapter::getReducedContinuousFunctionCode() const
146 {
147 return ReducedNodes;
148 }
149
150 int MeshAdapter::getFunctionCode() const
151 {
152 return Elements;
153 }
154 int MeshAdapter::getReducedFunctionCode() const
155 {
156 return ReducedElements;
157 }
158
159 int MeshAdapter::getFunctionOnBoundaryCode() const
160 {
161 return FaceElements;
162 }
163 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
164 {
165 return ReducedFaceElements;
166 }
167
168 int MeshAdapter::getFunctionOnContactZeroCode() const
169 {
170 return ContactElementsZero;
171 }
172 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
173 {
174 return ReducedContactElementsZero;
175 }
176
177 int MeshAdapter::getFunctionOnContactOneCode() const
178 {
179 return ContactElementsOne;
180 }
181 int MeshAdapter::getReducedFunctionOnContactOneCode() const
182 {
183 return ReducedContactElementsOne;
184 }
185
186 int MeshAdapter::getSolutionCode() const
187 {
188 return DegreesOfFreedom;
189 }
190
191 int MeshAdapter::getReducedSolutionCode() const
192 {
193 return ReducedDegreesOfFreedom;
194 }
195
196 int MeshAdapter::getDiracDeltaFunctionCode() const
197 {
198 return Points;
199 }
200
201 //
202 // return the spatial dimension of the Mesh:
203 //
204 int MeshAdapter::getDim() const
205 {
206 int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
207 checkFinleyError();
208 return numDim;
209 }
210
211 //
212 // return the number of data points per sample and the number of samples
213 // needed to represent data on a parts of the mesh.
214 //
215 pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
216 {
217 int numDataPointsPerSample=0;
218 int numSamples=0;
219 Finley_Mesh* mesh=m_finleyMesh.get();
220 switch (functionSpaceCode) {
221 case(Nodes):
222 numDataPointsPerSample=1;
223 if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
224 break;
225 case(ReducedNodes):
226 /* TODO: add ReducedNodes */
227 numDataPointsPerSample=1;
228 if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
229 throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");
230 break;
231 case(Elements):
232 if (mesh->Elements!=NULL) {
233 numSamples=mesh->Elements->numElements;
234 numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
235 }
236 break;
237 case(ReducedElements):
238 if (mesh->Elements!=NULL) {
239 numSamples=mesh->Elements->numElements;
240 numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
241 }
242 break;
243 case(FaceElements):
244 if (mesh->FaceElements!=NULL) {
245 numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
246 numSamples=mesh->FaceElements->numElements;
247 }
248 break;
249 case(ReducedFaceElements):
250 if (mesh->FaceElements!=NULL) {
251 numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
252 numSamples=mesh->FaceElements->numElements;
253 }
254 break;
255 case(Points):
256 if (mesh->Points!=NULL) {
257 numDataPointsPerSample=1;
258 numSamples=mesh->Points->numElements;
259 }
260 break;
261 case(ContactElementsZero):
262 if (mesh->ContactElements!=NULL) {
263 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
264 numSamples=mesh->ContactElements->numElements;
265 }
266 break;
267 case(ReducedContactElementsZero):
268 if (mesh->ContactElements!=NULL) {
269 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
270 numSamples=mesh->ContactElements->numElements;
271 }
272 break;
273 case(ContactElementsOne):
274 if (mesh->ContactElements!=NULL) {
275 numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
276 numSamples=mesh->ContactElements->numElements;
277 }
278 break;
279 case(ReducedContactElementsOne):
280 if (mesh->ContactElements!=NULL) {
281 numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
282 numSamples=mesh->ContactElements->numElements;
283 }
284 break;
285 case(DegreesOfFreedom):
286 if (mesh->Nodes!=NULL) {
287 numDataPointsPerSample=1;
288 #ifndef PASO_MPI
289 numSamples=mesh->Nodes->numDegreesOfFreedom;
290 #else
291 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
292 #endif
293 }
294 break;
295 case(ReducedDegreesOfFreedom):
296 if (mesh->Nodes!=NULL) {
297 numDataPointsPerSample=1;
298 #ifndef PASO_MPI
299 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
300 #else
301 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
302 #endif
303 }
304 break;
305 default:
306 stringstream temp;
307 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
308 throw FinleyAdapterException(temp.str());
309 break;
310 }
311 return pair<int,int>(numDataPointsPerSample,numSamples);
312 }
313
314 //
315 // adds linear PDE of second order into a given stiffness matrix and right hand side:
316 //
317 void MeshAdapter::addPDEToSystem(
318 SystemMatrixAdapter& mat, escript::Data& rhs,
319 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
320 const escript::Data& d, const escript::Data& y,
321 const escript::Data& d_contact,const escript::Data& y_contact) const
322 {
323 escriptDataC _rhs=rhs.getDataC();
324 escriptDataC _A =A.getDataC();
325 escriptDataC _B=B.getDataC();
326 escriptDataC _C=C.getDataC();
327 escriptDataC _D=D.getDataC();
328 escriptDataC _X=X.getDataC();
329 escriptDataC _Y=Y.getDataC();
330 escriptDataC _d=d.getDataC();
331 escriptDataC _y=y.getDataC();
332 escriptDataC _d_contact=d_contact.getDataC();
333 escriptDataC _y_contact=y_contact.getDataC();
334
335 Finley_Mesh* mesh=m_finleyMesh.get();
336
337 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
338 checkFinleyError();
339
340 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
341 checkFinleyError();
342
343 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
344 checkFinleyError();
345 }
346
347 //
348 // adds linear PDE of second order into the right hand side only
349 //
350 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
351 {
352 Finley_Mesh* mesh=m_finleyMesh.get();
353
354 escriptDataC _rhs=rhs.getDataC();
355 escriptDataC _X=X.getDataC();
356 escriptDataC _Y=Y.getDataC();
357 escriptDataC _y=y.getDataC();
358 escriptDataC _y_contact=y_contact.getDataC();
359
360 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
361 checkFinleyError();
362
363 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
364 checkFinleyError();
365
366 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
367 checkFinleyError();
368 }
369
370 //
371 // interpolates data between different function spaces:
372 //
373 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
374 {
375 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
376 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
377 if (inDomain!=*this)
378 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
379 if (targetDomain!=*this)
380 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
381
382 Finley_Mesh* mesh=m_finleyMesh.get();
383 escriptDataC _target=target.getDataC();
384 escriptDataC _in=in.getDataC();
385 switch(in.getFunctionSpace().getTypeCode()) {
386 case(Nodes):
387 switch(target.getFunctionSpace().getTypeCode()) {
388 case(Nodes):
389 case(ReducedNodes):
390 case(DegreesOfFreedom):
391 case(ReducedDegreesOfFreedom):
392 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
393 break;
394 case(Elements):
395 case(ReducedElements):
396 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
397 break;
398 case(FaceElements):
399 case(ReducedFaceElements):
400 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
401 break;
402 case(Points):
403 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
404 break;
405 case(ContactElementsZero):
406 case(ReducedContactElementsZero):
407 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
408 break;
409 case(ContactElementsOne):
410 case(ReducedContactElementsOne):
411 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
412 break;
413 default:
414 stringstream temp;
415 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
416 throw FinleyAdapterException(temp.str());
417 break;
418 }
419 break;
420 case(ReducedNodes):
421 switch(target.getFunctionSpace().getTypeCode()) {
422 case(Nodes):
423 case(ReducedNodes):
424 case(DegreesOfFreedom):
425 case(ReducedDegreesOfFreedom):
426 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
427 break;
428 case(Elements):
429 case(ReducedElements):
430 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
431 break;
432 case(FaceElements):
433 case(ReducedFaceElements):
434 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
435 break;
436 case(Points):
437 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
438 break;
439 case(ContactElementsZero):
440 case(ReducedContactElementsZero):
441 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
442 break;
443 case(ContactElementsOne):
444 case(ReducedContactElementsOne):
445 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
446 break;
447 default:
448 stringstream temp;
449 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
450 throw FinleyAdapterException(temp.str());
451 break;
452 }
453 break;
454 case(Elements):
455 if (target.getFunctionSpace().getTypeCode()==Elements) {
456 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
457 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
458 Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
459 } else {
460 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
461 }
462 break;
463 case(ReducedElements):
464 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
465 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
466 } else {
467 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
468 }
469 break;
470 case(FaceElements):
471 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
472 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
473 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
474 Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
475 } else {
476 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
477 }
478 break;
479 case(ReducedFaceElements):
480 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
481 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
482 } else {
483 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
484 }
485 break;
486 case(Points):
487 if (target.getFunctionSpace().getTypeCode()==Points) {
488 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
489 } else {
490 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
491 }
492 break;
493 case(ContactElementsZero):
494 case(ContactElementsOne):
495 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
496 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
497 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
498 Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
499 } else {
500 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
501 break;
502 }
503 break;
504 case(ReducedContactElementsZero):
505 case(ReducedContactElementsOne):
506 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
507 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
508 } else {
509 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
510 break;
511 }
512 break;
513 case(DegreesOfFreedom):
514 switch(target.getFunctionSpace().getTypeCode()) {
515 case(ReducedDegreesOfFreedom):
516 case(DegreesOfFreedom):
517 case(Nodes):
518 case(ReducedNodes):
519 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
520 break;
521 #ifndef PASO_MPI
522 case(Elements):
523 case(ReducedElements):
524 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
525 break;
526 case(FaceElements):
527 case(ReducedFaceElements):
528 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
529 break;
530 case(Points):
531 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
532 break;
533 case(ContactElementsZero):
534 case(ContactElementsOne):
535 case(ReducedContactElementsZero):
536 case(ReducedContactElementsOne):
537 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
538 break;
539 #else
540 /* need to copy Degrees of freedom data to nodal data so that the external values are available */
541 case(Elements):
542 case(ReducedElements):
543 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
544 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);
545 break;
546 case(FaceElements):
547 case(ReducedFaceElements):
548 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
549 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);
550 break;
551 case(Points):
552 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
553 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);
554 break;
555 case(ContactElementsZero):
556 case(ContactElementsOne):
557 case(ReducedContactElementsZero):
558 case(ReducedContactElementsOne):
559 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
560 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);
561 break;
562 #endif
563 default:
564 stringstream temp;
565 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
566 throw FinleyAdapterException(temp.str());
567 break;
568 }
569 break;
570 case(ReducedDegreesOfFreedom):
571 switch(target.getFunctionSpace().getTypeCode()) {
572 case(Nodes):
573 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
574 break;
575 case(ReducedNodes):
576 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
577 break;
578 case(DegreesOfFreedom):
579 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
580 break;
581 case(ReducedDegreesOfFreedom):
582 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
583 break;
584 case(Elements):
585 case(ReducedElements):
586 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
587 break;
588 case(FaceElements):
589 case(ReducedFaceElements):
590 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
591 break;
592 case(Points):
593 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
594 break;
595 case(ContactElementsZero):
596 case(ContactElementsOne):
597 case(ReducedContactElementsZero):
598 case(ReducedContactElementsOne):
599 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
600 break;
601 default:
602 stringstream temp;
603 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
604 throw FinleyAdapterException(temp.str());
605 break;
606 }
607 break;
608 default:
609 stringstream temp;
610 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
611 throw FinleyAdapterException(temp.str());
612 break;
613 }
614 checkFinleyError();
615 }
616
617 //
618 // copies the locations of sample points into x:
619 //
620 void MeshAdapter::setToX(escript::Data& arg) const
621 {
622 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
623 if (argDomain!=*this)
624 throw FinleyAdapterException("Error - Illegal domain of data point locations");
625 Finley_Mesh* mesh=m_finleyMesh.get();
626 // in case of values node coordinates we can do the job directly:
627 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
628 escriptDataC _arg=arg.getDataC();
629 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
630 } else {
631 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
632 escriptDataC _tmp_data=tmp_data.getDataC();
633 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
634 // this is then interpolated onto arg:
635 interpolateOnDomain(arg,tmp_data);
636 }
637 checkFinleyError();
638 }
639
640 //
641 // return the normal vectors at the location of data points as a Data object:
642 //
643 void MeshAdapter::setToNormal(escript::Data& normal) const
644 {
645 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
646 if (normalDomain!=*this)
647 throw FinleyAdapterException("Error - Illegal domain of normal locations");
648 Finley_Mesh* mesh=m_finleyMesh.get();
649 escriptDataC _normal=normal.getDataC();
650 switch(normal.getFunctionSpace().getTypeCode()) {
651 case(Nodes):
652 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
653 break;
654 case(ReducedNodes):
655 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
656 break;
657 case(Elements):
658 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
659 break;
660 case(ReducedElements):
661 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
662 break;
663 case (FaceElements):
664 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
665 break;
666 case (ReducedFaceElements):
667 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
668 break;
669 case(Points):
670 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
671 break;
672 case (ContactElementsOne):
673 case (ContactElementsZero):
674 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
675 break;
676 case (ReducedContactElementsOne):
677 case (ReducedContactElementsZero):
678 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
679 break;
680 case(DegreesOfFreedom):
681 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
682 break;
683 case(ReducedDegreesOfFreedom):
684 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
685 break;
686 default:
687 stringstream temp;
688 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
689 throw FinleyAdapterException(temp.str());
690 break;
691 }
692 checkFinleyError();
693 }
694
695 //
696 // interpolates data to other domain:
697 //
698 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
699 {
700 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
701 if (targetDomain!=*this)
702 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
703
704 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
705 }
706
707 //
708 // calculates the integral of a function defined of arg:
709 //
710 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
711 {
712 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
713 if (argDomain!=*this)
714 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
715
716 Finley_Mesh* mesh=m_finleyMesh.get();
717 escriptDataC _arg=arg.getDataC();
718 switch(arg.getFunctionSpace().getTypeCode()) {
719 case(Nodes):
720 /* TODO */
721 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
722 break;
723 case(ReducedNodes):
724 /* TODO */
725 throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");
726 break;
727 case(Elements):
728 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
729 break;
730 case(ReducedElements):
731 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
732 break;
733 case(FaceElements):
734 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
735 break;
736 case(ReducedFaceElements):
737 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
738 break;
739 case(Points):
740 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
741 break;
742 case(ContactElementsZero):
743 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
744 break;
745 case(ReducedContactElementsZero):
746 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
747 break;
748 case(ContactElementsOne):
749 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
750 break;
751 case(ReducedContactElementsOne):
752 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
753 break;
754 case(DegreesOfFreedom):
755 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
756 break;
757 case(ReducedDegreesOfFreedom):
758 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
759 break;
760 default:
761 stringstream temp;
762 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
763 throw FinleyAdapterException(temp.str());
764 break;
765 }
766 checkFinleyError();
767 }
768
769 //
770 // calculates the gradient of arg:
771 //
772 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
773 {
774 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
775 if (argDomain!=*this)
776 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
777 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
778 if (gradDomain!=*this)
779 throw FinleyAdapterException("Error - Illegal domain of gradient");
780
781 Finley_Mesh* mesh=m_finleyMesh.get();
782 escriptDataC _grad=grad.getDataC();
783 escriptDataC nodeDataC;
784 #ifdef PASO_MPI
785 escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
786 if( arg.getFunctionSpace().getTypeCode() != Nodes )
787 {
788 Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
789 nodeDataC = nodeTemp.getDataC();
790 }
791 else
792 nodeDataC = arg.getDataC();
793 #else
794 nodeDataC = arg.getDataC();
795 #endif
796 switch(grad.getFunctionSpace().getTypeCode()) {
797 case(Nodes):
798 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
799 break;
800 case(ReducedNodes):
801 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
802 break;
803 case(Elements):
804 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
805 break;
806 case(ReducedElements):
807 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
808 break;
809 case(FaceElements):
810 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
811 break;
812 case(ReducedFaceElements):
813 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
814 break;
815 case(Points):
816 throw FinleyAdapterException("Error - Gradient at points is not supported.");
817 break;
818 case(ContactElementsZero):
819 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
820 break;
821 case(ReducedContactElementsZero):
822 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
823 break;
824 case(ContactElementsOne):
825 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
826 break;
827 case(ReducedContactElementsOne):
828 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
829 break;
830 case(DegreesOfFreedom):
831 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
832 break;
833 case(ReducedDegreesOfFreedom):
834 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
835 break;
836 default:
837 stringstream temp;
838 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
839 throw FinleyAdapterException(temp.str());
840 break;
841 }
842 checkFinleyError();
843 }
844
845 //
846 // returns the size of elements:
847 //
848 void MeshAdapter::setToSize(escript::Data& size) const
849 {
850 Finley_Mesh* mesh=m_finleyMesh.get();
851 escriptDataC tmp=size.getDataC();
852 switch(size.getFunctionSpace().getTypeCode()) {
853 case(Nodes):
854 throw FinleyAdapterException("Error - Size of nodes is not supported.");
855 break;
856 case(ReducedNodes):
857 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
858 break;
859 case(Elements):
860 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
861 break;
862 case(ReducedElements):
863 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
864 break;
865 case(FaceElements):
866 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
867 break;
868 case(ReducedFaceElements):
869 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
870 break;
871 case(Points):
872 throw FinleyAdapterException("Error - Size of point elements is not supported.");
873 break;
874 case(ContactElementsZero):
875 case(ContactElementsOne):
876 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
877 break;
878 case(ReducedContactElementsZero):
879 case(ReducedContactElementsOne):
880 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
881 break;
882 case(DegreesOfFreedom):
883 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
884 break;
885 case(ReducedDegreesOfFreedom):
886 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
887 break;
888 default:
889 stringstream temp;
890 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
891 throw FinleyAdapterException(temp.str());
892 break;
893 }
894 checkFinleyError();
895 }
896
897 // sets the location of nodes:
898 void MeshAdapter::setNewX(const escript::Data& new_x)
899 {
900 Finley_Mesh* mesh=m_finleyMesh.get();
901 escriptDataC tmp;
902 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
903 if (newDomain!=*this)
904 throw FinleyAdapterException("Error - Illegal domain of new point locations");
905 tmp = new_x.getDataC();
906 Finley_Mesh_setCoordinates(mesh,&tmp);
907 checkFinleyError();
908 }
909
910 // saves a data array in openDX format:
911 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
912 {
913 int MAX_namelength=256;
914 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
915 /* win32 refactor */
916 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
917 for(int i=0;i<num_data;i++)
918 {
919 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
920 }
921
922 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
923 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
924 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
925
926 boost::python::list keys=arg.keys();
927 for (int i=0;i<num_data;++i) {
928 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
929 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
930 throw FinleyAdapterException("Error in saveDX: Data must be defined on same Domain");
931 data[i]=d.getDataC();
932 ptr_data[i]=&(data[i]);
933 std::string n=boost::python::extract<std::string>(keys[i]);
934 c_names[i]=&(names[i][0]);
935 if (n.length()>MAX_namelength-1) {
936 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
937 c_names[i][MAX_namelength-1]='\0';
938 } else {
939 strcpy(c_names[i],n.c_str());
940 }
941 }
942 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
943 checkFinleyError();
944
945 /* win32 refactor */
946 TMPMEMFREE(c_names);
947 TMPMEMFREE(data);
948 TMPMEMFREE(ptr_data);
949 for(int i=0;i<num_data;i++)
950 {
951 TMPMEMFREE(names[i]);
952 }
953 TMPMEMFREE(names);
954
955 return;
956 }
957
958 // saves a data array in openVTK format:
959 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
960 {
961 int MAX_namelength=256;
962 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
963 /* win32 refactor */
964 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
965 for(int i=0;i<num_data;i++)
966 {
967 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
968 }
969
970 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
971 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
972 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
973
974 boost::python::list keys=arg.keys();
975 for (int i=0;i<num_data;++i) {
976 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
977 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
978 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
979 data[i]=d.getDataC();
980 ptr_data[i]=&(data[i]);
981 std::string n=boost::python::extract<std::string>(keys[i]);
982 c_names[i]=&(names[i][0]);
983 if (n.length()>MAX_namelength-1) {
984 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
985 c_names[i][MAX_namelength-1]='\0';
986 } else {
987 strcpy(c_names[i],n.c_str());
988 }
989 }
990 #ifndef PASO_MPI
991 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
992 #else
993 Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
994 #endif
995
996 checkFinleyError();
997 /* win32 refactor */
998 TMPMEMFREE(c_names);
999 TMPMEMFREE(data);
1000 TMPMEMFREE(ptr_data);
1001 for(int i=0;i<num_data;i++)
1002 {
1003 TMPMEMFREE(names[i]);
1004 }
1005 TMPMEMFREE(names);
1006
1007 return;
1008 }
1009
1010
1011 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1012 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1013 const int row_blocksize,
1014 const escript::FunctionSpace& row_functionspace,
1015 const int column_blocksize,
1016 const escript::FunctionSpace& column_functionspace,
1017 const int type) const
1018 {
1019 int reduceRowOrder=0;
1020 int reduceColOrder=0;
1021 // is the domain right?
1022 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1023 if (row_domain!=*this)
1024 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1025 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1026 if (col_domain!=*this)
1027 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1028 // is the function space type right
1029 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1030 reduceRowOrder=0;
1031 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1032 reduceRowOrder=1;
1033 } else {
1034 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1035 }
1036 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1037 reduceColOrder=0;
1038 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1039 reduceColOrder=1;
1040 } else {
1041 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1042 }
1043 // generate matrix:
1044
1045 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1046 checkFinleyError();
1047 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1048 checkPasoError();
1049 Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
1050 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1051 }
1052
1053 //
1054 // vtkObject MeshAdapter::createVtkObject() const
1055 // TODO:
1056 //
1057
1058 //
1059 // returns true if data at the atom_type is considered as being cell centered:
1060 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1061 {
1062 switch(functionSpaceCode) {
1063 case(Nodes):
1064 case(DegreesOfFreedom):
1065 case(ReducedDegreesOfFreedom):
1066 return false;
1067 break;
1068 case(Elements):
1069 case(FaceElements):
1070 case(Points):
1071 case(ContactElementsZero):
1072 case(ContactElementsOne):
1073 case(ReducedElements):
1074 case(ReducedFaceElements):
1075 case(ReducedContactElementsZero):
1076 case(ReducedContactElementsOne):
1077 return true;
1078 break;
1079 default:
1080 stringstream temp;
1081 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1082 throw FinleyAdapterException(temp.str());
1083 break;
1084 }
1085 checkFinleyError();
1086 return false;
1087 }
1088
1089 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1090 {
1091 switch(functionSpaceType_source) {
1092 case(Nodes):
1093 switch(functionSpaceType_target) {
1094 case(Nodes):
1095 case(ReducedNodes):
1096 case(ReducedDegreesOfFreedom):
1097 case(DegreesOfFreedom):
1098 case(Elements):
1099 case(ReducedElements):
1100 case(FaceElements):
1101 case(ReducedFaceElements):
1102 case(Points):
1103 case(ContactElementsZero):
1104 case(ReducedContactElementsZero):
1105 case(ContactElementsOne):
1106 case(ReducedContactElementsOne):
1107 return true;
1108 default:
1109 stringstream temp;
1110 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1111 throw FinleyAdapterException(temp.str());
1112 }
1113 break;
1114 case(ReducedNodes):
1115 switch(functionSpaceType_target) {
1116 case(ReducedNodes):
1117 case(ReducedDegreesOfFreedom):
1118 case(Elements):
1119 case(ReducedElements):
1120 case(FaceElements):
1121 case(ReducedFaceElements):
1122 case(Points):
1123 case(ContactElementsZero):
1124 case(ReducedContactElementsZero):
1125 case(ContactElementsOne):
1126 case(ReducedContactElementsOne):
1127 return true;
1128 case(Nodes):
1129 case(DegreesOfFreedom):
1130 return false;
1131 default:
1132 stringstream temp;
1133 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1134 throw FinleyAdapterException(temp.str());
1135 }
1136 break;
1137 case(Elements):
1138 if (functionSpaceType_target==Elements) {
1139 return true;
1140 } else if (functionSpaceType_target==ReducedElements) {
1141 return true;
1142 } else {
1143 return false;
1144 }
1145 case(ReducedElements):
1146 if (functionSpaceType_target==ReducedElements) {
1147 return true;
1148 } else {
1149 return false;
1150 }
1151 case(FaceElements):
1152 if (functionSpaceType_target==FaceElements) {
1153 return true;
1154 } else if (functionSpaceType_target==ReducedFaceElements) {
1155 return true;
1156 } else {
1157 return false;
1158 }
1159 case(ReducedFaceElements):
1160 if (functionSpaceType_target==ReducedFaceElements) {
1161 return true;
1162 } else {
1163 return false;
1164 }
1165 case(Points):
1166 if (functionSpaceType_target==Points) {
1167 return true;
1168 } else {
1169 return false;
1170 }
1171 case(ContactElementsZero):
1172 case(ContactElementsOne):
1173 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1174 return true;
1175 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1176 return true;
1177 } else {
1178 return false;
1179 }
1180 case(ReducedContactElementsZero):
1181 case(ReducedContactElementsOne):
1182 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1183 return true;
1184 } else {
1185 return false;
1186 }
1187 case(DegreesOfFreedom):
1188 switch(functionSpaceType_target) {
1189 case(ReducedDegreesOfFreedom):
1190 case(DegreesOfFreedom):
1191 case(Nodes):
1192 case(ReducedNodes):
1193 case(Elements):
1194 case(ReducedElements):
1195 case(Points):
1196 case(FaceElements):
1197 case(ReducedFaceElements):
1198 case(ContactElementsZero):
1199 case(ReducedContactElementsZero):
1200 case(ContactElementsOne):
1201 case(ReducedContactElementsOne):
1202 return true;
1203 default:
1204 stringstream temp;
1205 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1206 throw FinleyAdapterException(temp.str());
1207 }
1208 break;
1209 case(ReducedDegreesOfFreedom):
1210 switch(functionSpaceType_target) {
1211 case(ReducedDegreesOfFreedom):
1212 case(ReducedNodes):
1213 case(Elements):
1214 case(ReducedElements):
1215 case(FaceElements):
1216 case(ReducedFaceElements):
1217 case(Points):
1218 case(ContactElementsZero):
1219 case(ReducedContactElementsZero):
1220 case(ContactElementsOne):
1221 case(ReducedContactElementsOne):
1222 return true;
1223 case(Nodes):
1224 case(DegreesOfFreedom):
1225 return false;
1226 default:
1227 stringstream temp;
1228 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1229 throw FinleyAdapterException(temp.str());
1230 }
1231 break;
1232 default:
1233 stringstream temp;
1234 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1235 throw FinleyAdapterException(temp.str());
1236 break;
1237 }
1238 checkFinleyError();
1239 return false;
1240 }
1241
1242 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1243 {
1244 return false;
1245 }
1246
1247 bool MeshAdapter::operator==(const AbstractDomain& other) const
1248 {
1249 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1250 if (temp!=0) {
1251 return (m_finleyMesh==temp->m_finleyMesh);
1252 } else {
1253 return false;
1254 }
1255 }
1256
1257 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1258 {
1259 return !(operator==(other));
1260 }
1261
1262 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1263 {
1264 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1265 checkPasoError();
1266 return out;
1267 }
1268
1269 escript::Data MeshAdapter::getX() const
1270 {
1271 return continuousFunction(asAbstractContinuousDomain()).getX();
1272 }
1273
1274 escript::Data MeshAdapter::getNormal() const
1275 {
1276 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1277 }
1278
1279 escript::Data MeshAdapter::getSize() const
1280 {
1281 return function(asAbstractContinuousDomain()).getSize();
1282 }
1283
1284 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1285 {
1286 int *out=0,i;
1287 Finley_Mesh* mesh=m_finleyMesh.get();
1288 switch (functionSpaceType) {
1289 case(Nodes):
1290 if (mesh->Nodes!=NULL) {
1291 out=mesh->Nodes->Id;
1292 break;
1293 }
1294 case(ReducedNodes):
1295 throw FinleyAdapterException("Error - ReducedNodes not supported yet.");
1296 break;
1297 case(Elements):
1298 out=mesh->Elements->Id;
1299 break;
1300 case(ReducedElements):
1301 out=mesh->Elements->Id;
1302 break;
1303 case(FaceElements):
1304 out=mesh->FaceElements->Id;
1305 break;
1306 case(ReducedFaceElements):
1307 out=mesh->FaceElements->Id;
1308 break;
1309 case(Points):
1310 out=mesh->Points->Id;
1311 break;
1312 case(ContactElementsZero):
1313 out=mesh->ContactElements->Id;
1314 break;
1315 case(ReducedContactElementsZero):
1316 out=mesh->ContactElements->Id;
1317 break;
1318 case(ContactElementsOne):
1319 out=mesh->ContactElements->Id;
1320 break;
1321 case(ReducedContactElementsOne):
1322 out=mesh->ContactElements->Id;
1323 break;
1324 case(DegreesOfFreedom):
1325 out=mesh->Nodes->degreeOfFreedomId;
1326 break;
1327 case(ReducedDegreesOfFreedom):
1328 out=mesh->Nodes->reducedDegreeOfFreedomId;
1329 break;
1330 default:
1331 stringstream temp;
1332 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1333 throw FinleyAdapterException(temp.str());
1334 break;
1335 }
1336 return out;
1337 }
1338 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1339 {
1340 int out=0;
1341 Finley_Mesh* mesh=m_finleyMesh.get();
1342 switch (functionSpaceType) {
1343 case(Nodes):
1344 out=mesh->Nodes->Tag[sampleNo];
1345 break;
1346 case(ReducedNodes):
1347 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1348 break;
1349 case(Elements):
1350 out=mesh->Elements->Tag[sampleNo];
1351 break;
1352 case(ReducedElements):
1353 out=mesh->Elements->Tag[sampleNo];
1354 break;
1355 case(FaceElements):
1356 out=mesh->FaceElements->Tag[sampleNo];
1357 break;
1358 case(ReducedFaceElements):
1359 out=mesh->FaceElements->Tag[sampleNo];
1360 break;
1361 case(Points):
1362 out=mesh->Points->Tag[sampleNo];
1363 break;
1364 case(ContactElementsZero):
1365 out=mesh->ContactElements->Tag[sampleNo];
1366 break;
1367 case(ReducedContactElementsZero):
1368 out=mesh->ContactElements->Tag[sampleNo];
1369 break;
1370 case(ContactElementsOne):
1371 out=mesh->ContactElements->Tag[sampleNo];
1372 break;
1373 case(ReducedContactElementsOne):
1374 out=mesh->ContactElements->Tag[sampleNo];
1375 break;
1376 case(DegreesOfFreedom):
1377 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1378 break;
1379 case(ReducedDegreesOfFreedom):
1380 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1381 break;
1382 default:
1383 stringstream temp;
1384 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1385 throw FinleyAdapterException(temp.str());
1386 break;
1387 }
1388 return out;
1389 }
1390
1391
1392 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1393 {
1394 Finley_Mesh* mesh=m_finleyMesh.get();
1395 escriptDataC tmp=mask.getDataC();
1396 switch(functionSpaceType) {
1397 case(Nodes):
1398 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1399 break;
1400 case(ReducedNodes):
1401 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1402 break;
1403 case(DegreesOfFreedom):
1404 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1405 break;
1406 case(ReducedDegreesOfFreedom):
1407 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1408 break;
1409 case(Elements):
1410 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1411 break;
1412 case(ReducedElements):
1413 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1414 break;
1415 case(FaceElements):
1416 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1417 break;
1418 case(ReducedFaceElements):
1419 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1420 break;
1421 case(Points):
1422 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1423 break;
1424 case(ContactElementsZero):
1425 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1426 break;
1427 case(ReducedContactElementsZero):
1428 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1429 break;
1430 case(ContactElementsOne):
1431 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1432 break;
1433 case(ReducedContactElementsOne):
1434 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1435 break;
1436 default:
1437 stringstream temp;
1438 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1439 throw FinleyAdapterException(temp.str());
1440 }
1441 checkFinleyError();
1442 return;
1443 }
1444
1445 void MeshAdapter::setTagMap(const std::string& name, int tag)
1446 {
1447 Finley_Mesh* mesh=m_finleyMesh.get();
1448 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1449 checkPasoError();
1450 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1451 }
1452
1453 int MeshAdapter::getTag(const std::string& name) const
1454 {
1455 Finley_Mesh* mesh=m_finleyMesh.get();
1456 int tag=0;
1457 tag=Finley_Mesh_getTag(mesh, name.c_str());
1458 checkPasoError();
1459 // throwStandardException("MeshAdapter::getTag is not implemented.");
1460 return tag;
1461 }
1462
1463 bool MeshAdapter::isValidTagName(const std::string& name) const
1464 {
1465 Finley_Mesh* mesh=m_finleyMesh.get();
1466 return Finley_Mesh_isValidTagName(mesh,name.c_str());
1467 }
1468
1469 std::string MeshAdapter::showTagNames() const
1470 {
1471 stringstream temp;
1472 Finley_Mesh* mesh=m_finleyMesh.get();
1473 Finley_TagMap* tag_map=mesh->TagMap;
1474 while (tag_map) {
1475 temp << tag_map->name;
1476 tag_map=tag_map->next;
1477 if (tag_map) temp << ", ";
1478 }
1479 return string(temp.str());
1480 }
1481
1482 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26