/[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 1204 - (show annotations)
Sat Jun 23 11:43:12 2007 UTC (11 years, 11 months ago) by gross
File size: 56656 byte(s)
a frame for an improved lumping procedure added.
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 void MeshAdapter::addPDEToLumpedSystem(
348 escript::Data& mat,
349 const escript::Data& D,
350 const escript::Data& d) const
351 {
352 escriptDataC _mat=mat.getDataC();
353 escriptDataC _D=D.getDataC();
354 escriptDataC _d=d.getDataC();
355
356 Finley_Mesh* mesh=m_finleyMesh.get();
357
358 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
359 Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
360
361 checkFinleyError();
362 }
363
364
365 //
366 // adds linear PDE of second order into the right hand side only
367 //
368 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
369 {
370 Finley_Mesh* mesh=m_finleyMesh.get();
371
372 escriptDataC _rhs=rhs.getDataC();
373 escriptDataC _X=X.getDataC();
374 escriptDataC _Y=Y.getDataC();
375 escriptDataC _y=y.getDataC();
376 escriptDataC _y_contact=y_contact.getDataC();
377
378 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
379 checkFinleyError();
380
381 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
382 checkFinleyError();
383
384 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
385 checkFinleyError();
386 }
387
388 //
389 // interpolates data between different function spaces:
390 //
391 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
392 {
393 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
394 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
395 if (inDomain!=*this)
396 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
397 if (targetDomain!=*this)
398 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
399
400 Finley_Mesh* mesh=m_finleyMesh.get();
401 escriptDataC _target=target.getDataC();
402 escriptDataC _in=in.getDataC();
403 switch(in.getFunctionSpace().getTypeCode()) {
404 case(Nodes):
405 switch(target.getFunctionSpace().getTypeCode()) {
406 case(Nodes):
407 case(ReducedNodes):
408 case(DegreesOfFreedom):
409 case(ReducedDegreesOfFreedom):
410 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
411 break;
412 case(Elements):
413 case(ReducedElements):
414 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
415 break;
416 case(FaceElements):
417 case(ReducedFaceElements):
418 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
419 break;
420 case(Points):
421 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
422 break;
423 case(ContactElementsZero):
424 case(ReducedContactElementsZero):
425 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
426 break;
427 case(ContactElementsOne):
428 case(ReducedContactElementsOne):
429 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
430 break;
431 default:
432 stringstream temp;
433 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
434 throw FinleyAdapterException(temp.str());
435 break;
436 }
437 break;
438 case(ReducedNodes):
439 switch(target.getFunctionSpace().getTypeCode()) {
440 case(Nodes):
441 case(ReducedNodes):
442 case(DegreesOfFreedom):
443 case(ReducedDegreesOfFreedom):
444 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
445 break;
446 case(Elements):
447 case(ReducedElements):
448 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
449 break;
450 case(FaceElements):
451 case(ReducedFaceElements):
452 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
453 break;
454 case(Points):
455 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
456 break;
457 case(ContactElementsZero):
458 case(ReducedContactElementsZero):
459 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
460 break;
461 case(ContactElementsOne):
462 case(ReducedContactElementsOne):
463 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
464 break;
465 default:
466 stringstream temp;
467 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
468 throw FinleyAdapterException(temp.str());
469 break;
470 }
471 break;
472 case(Elements):
473 if (target.getFunctionSpace().getTypeCode()==Elements) {
474 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
475 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
476 Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
477 } else {
478 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
479 }
480 break;
481 case(ReducedElements):
482 if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
483 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
484 } else {
485 throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
486 }
487 break;
488 case(FaceElements):
489 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
490 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
491 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
492 Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
493 } else {
494 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
495 }
496 break;
497 case(ReducedFaceElements):
498 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
499 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
500 } else {
501 throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
502 }
503 break;
504 case(Points):
505 if (target.getFunctionSpace().getTypeCode()==Points) {
506 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
507 } else {
508 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
509 }
510 break;
511 case(ContactElementsZero):
512 case(ContactElementsOne):
513 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
514 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
515 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
516 Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
517 } else {
518 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
519 break;
520 }
521 break;
522 case(ReducedContactElementsZero):
523 case(ReducedContactElementsOne):
524 if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
525 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
526 } else {
527 throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
528 break;
529 }
530 break;
531 case(DegreesOfFreedom):
532 switch(target.getFunctionSpace().getTypeCode()) {
533 case(ReducedDegreesOfFreedom):
534 case(DegreesOfFreedom):
535 case(Nodes):
536 case(ReducedNodes):
537 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
538 break;
539 #ifndef PASO_MPI
540 case(Elements):
541 case(ReducedElements):
542 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
543 break;
544 case(FaceElements):
545 case(ReducedFaceElements):
546 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
547 break;
548 case(Points):
549 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
550 break;
551 case(ContactElementsZero):
552 case(ContactElementsOne):
553 case(ReducedContactElementsZero):
554 case(ReducedContactElementsOne):
555 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
556 break;
557 #else
558 /* need to copy Degrees of freedom data to nodal data so that the external values are available */
559 case(Elements):
560 case(ReducedElements):
561 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
562 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);
563 break;
564 case(FaceElements):
565 case(ReducedFaceElements):
566 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
567 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);
568 break;
569 case(Points):
570 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
571 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);
572 break;
573 case(ContactElementsZero):
574 case(ContactElementsOne):
575 case(ReducedContactElementsZero):
576 case(ReducedContactElementsOne):
577 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
578 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);
579 break;
580 #endif
581 default:
582 stringstream temp;
583 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
584 throw FinleyAdapterException(temp.str());
585 break;
586 }
587 break;
588 case(ReducedDegreesOfFreedom):
589 switch(target.getFunctionSpace().getTypeCode()) {
590 case(Nodes):
591 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
592 break;
593 case(ReducedNodes):
594 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
595 break;
596 case(DegreesOfFreedom):
597 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
598 break;
599 case(ReducedDegreesOfFreedom):
600 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
601 break;
602 case(Elements):
603 case(ReducedElements):
604 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
605 break;
606 case(FaceElements):
607 case(ReducedFaceElements):
608 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
609 break;
610 case(Points):
611 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
612 break;
613 case(ContactElementsZero):
614 case(ContactElementsOne):
615 case(ReducedContactElementsZero):
616 case(ReducedContactElementsOne):
617 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
618 break;
619 default:
620 stringstream temp;
621 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
622 throw FinleyAdapterException(temp.str());
623 break;
624 }
625 break;
626 default:
627 stringstream temp;
628 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
629 throw FinleyAdapterException(temp.str());
630 break;
631 }
632 checkFinleyError();
633 }
634
635 //
636 // copies the locations of sample points into x:
637 //
638 void MeshAdapter::setToX(escript::Data& arg) const
639 {
640 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
641 if (argDomain!=*this)
642 throw FinleyAdapterException("Error - Illegal domain of data point locations");
643 Finley_Mesh* mesh=m_finleyMesh.get();
644 // in case of values node coordinates we can do the job directly:
645 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
646 escriptDataC _arg=arg.getDataC();
647 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
648 } else {
649 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
650 escriptDataC _tmp_data=tmp_data.getDataC();
651 Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
652 // this is then interpolated onto arg:
653 interpolateOnDomain(arg,tmp_data);
654 }
655 checkFinleyError();
656 }
657
658 //
659 // return the normal vectors at the location of data points as a Data object:
660 //
661 void MeshAdapter::setToNormal(escript::Data& normal) const
662 {
663 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
664 if (normalDomain!=*this)
665 throw FinleyAdapterException("Error - Illegal domain of normal locations");
666 Finley_Mesh* mesh=m_finleyMesh.get();
667 escriptDataC _normal=normal.getDataC();
668 switch(normal.getFunctionSpace().getTypeCode()) {
669 case(Nodes):
670 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
671 break;
672 case(ReducedNodes):
673 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
674 break;
675 case(Elements):
676 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
677 break;
678 case(ReducedElements):
679 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
680 break;
681 case (FaceElements):
682 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
683 break;
684 case (ReducedFaceElements):
685 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
686 break;
687 case(Points):
688 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
689 break;
690 case (ContactElementsOne):
691 case (ContactElementsZero):
692 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
693 break;
694 case (ReducedContactElementsOne):
695 case (ReducedContactElementsZero):
696 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
697 break;
698 case(DegreesOfFreedom):
699 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
700 break;
701 case(ReducedDegreesOfFreedom):
702 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
703 break;
704 default:
705 stringstream temp;
706 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
707 throw FinleyAdapterException(temp.str());
708 break;
709 }
710 checkFinleyError();
711 }
712
713 //
714 // interpolates data to other domain:
715 //
716 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
717 {
718 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
719 if (targetDomain!=*this)
720 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
721
722 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
723 }
724
725 //
726 // calculates the integral of a function defined of arg:
727 //
728 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
729 {
730 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
731 if (argDomain!=*this)
732 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
733
734 Finley_Mesh* mesh=m_finleyMesh.get();
735 escriptDataC _arg=arg.getDataC();
736 switch(arg.getFunctionSpace().getTypeCode()) {
737 case(Nodes):
738 /* TODO */
739 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
740 break;
741 case(ReducedNodes):
742 /* TODO */
743 throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");
744 break;
745 case(Elements):
746 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
747 break;
748 case(ReducedElements):
749 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
750 break;
751 case(FaceElements):
752 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
753 break;
754 case(ReducedFaceElements):
755 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
756 break;
757 case(Points):
758 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
759 break;
760 case(ContactElementsZero):
761 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
762 break;
763 case(ReducedContactElementsZero):
764 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
765 break;
766 case(ContactElementsOne):
767 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
768 break;
769 case(ReducedContactElementsOne):
770 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
771 break;
772 case(DegreesOfFreedom):
773 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
774 break;
775 case(ReducedDegreesOfFreedom):
776 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
777 break;
778 default:
779 stringstream temp;
780 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
781 throw FinleyAdapterException(temp.str());
782 break;
783 }
784 checkFinleyError();
785 }
786
787 //
788 // calculates the gradient of arg:
789 //
790 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
791 {
792 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
793 if (argDomain!=*this)
794 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
795 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
796 if (gradDomain!=*this)
797 throw FinleyAdapterException("Error - Illegal domain of gradient");
798
799 Finley_Mesh* mesh=m_finleyMesh.get();
800 escriptDataC _grad=grad.getDataC();
801 escriptDataC nodeDataC;
802 #ifdef PASO_MPI
803 escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
804 if( arg.getFunctionSpace().getTypeCode() != Nodes )
805 {
806 Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
807 nodeDataC = nodeTemp.getDataC();
808 }
809 else
810 nodeDataC = arg.getDataC();
811 #else
812 nodeDataC = arg.getDataC();
813 #endif
814 switch(grad.getFunctionSpace().getTypeCode()) {
815 case(Nodes):
816 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
817 break;
818 case(ReducedNodes):
819 throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
820 break;
821 case(Elements):
822 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
823 break;
824 case(ReducedElements):
825 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
826 break;
827 case(FaceElements):
828 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
829 break;
830 case(ReducedFaceElements):
831 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
832 break;
833 case(Points):
834 throw FinleyAdapterException("Error - Gradient at points is not supported.");
835 break;
836 case(ContactElementsZero):
837 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
838 break;
839 case(ReducedContactElementsZero):
840 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
841 break;
842 case(ContactElementsOne):
843 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
844 break;
845 case(ReducedContactElementsOne):
846 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
847 break;
848 case(DegreesOfFreedom):
849 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
850 break;
851 case(ReducedDegreesOfFreedom):
852 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
853 break;
854 default:
855 stringstream temp;
856 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
857 throw FinleyAdapterException(temp.str());
858 break;
859 }
860 checkFinleyError();
861 }
862
863 //
864 // returns the size of elements:
865 //
866 void MeshAdapter::setToSize(escript::Data& size) const
867 {
868 Finley_Mesh* mesh=m_finleyMesh.get();
869 escriptDataC tmp=size.getDataC();
870 switch(size.getFunctionSpace().getTypeCode()) {
871 case(Nodes):
872 throw FinleyAdapterException("Error - Size of nodes is not supported.");
873 break;
874 case(ReducedNodes):
875 throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
876 break;
877 case(Elements):
878 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
879 break;
880 case(ReducedElements):
881 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
882 break;
883 case(FaceElements):
884 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
885 break;
886 case(ReducedFaceElements):
887 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
888 break;
889 case(Points):
890 throw FinleyAdapterException("Error - Size of point elements is not supported.");
891 break;
892 case(ContactElementsZero):
893 case(ContactElementsOne):
894 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
895 break;
896 case(ReducedContactElementsZero):
897 case(ReducedContactElementsOne):
898 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
899 break;
900 case(DegreesOfFreedom):
901 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
902 break;
903 case(ReducedDegreesOfFreedom):
904 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
905 break;
906 default:
907 stringstream temp;
908 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
909 throw FinleyAdapterException(temp.str());
910 break;
911 }
912 checkFinleyError();
913 }
914
915 // sets the location of nodes:
916 void MeshAdapter::setNewX(const escript::Data& new_x)
917 {
918 Finley_Mesh* mesh=m_finleyMesh.get();
919 escriptDataC tmp;
920 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
921 if (newDomain!=*this)
922 throw FinleyAdapterException("Error - Illegal domain of new point locations");
923 tmp = new_x.getDataC();
924 Finley_Mesh_setCoordinates(mesh,&tmp);
925 checkFinleyError();
926 }
927
928 // saves a data array in openDX format:
929 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
930 {
931 int MAX_namelength=256;
932 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
933 /* win32 refactor */
934 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
935 for(int i=0;i<num_data;i++)
936 {
937 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
938 }
939
940 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
941 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
942 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
943
944 boost::python::list keys=arg.keys();
945 for (int i=0;i<num_data;++i) {
946 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
947 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
948 throw FinleyAdapterException("Error in saveDX: Data must be defined on same Domain");
949 data[i]=d.getDataC();
950 ptr_data[i]=&(data[i]);
951 std::string n=boost::python::extract<std::string>(keys[i]);
952 c_names[i]=&(names[i][0]);
953 if (n.length()>MAX_namelength-1) {
954 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
955 c_names[i][MAX_namelength-1]='\0';
956 } else {
957 strcpy(c_names[i],n.c_str());
958 }
959 }
960 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
961 checkFinleyError();
962
963 /* win32 refactor */
964 TMPMEMFREE(c_names);
965 TMPMEMFREE(data);
966 TMPMEMFREE(ptr_data);
967 for(int i=0;i<num_data;i++)
968 {
969 TMPMEMFREE(names[i]);
970 }
971 TMPMEMFREE(names);
972
973 return;
974 }
975
976 // saves a data array in openVTK format:
977 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
978 {
979 int MAX_namelength=256;
980 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
981 /* win32 refactor */
982 char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
983 for(int i=0;i<num_data;i++)
984 {
985 names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
986 }
987
988 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
989 escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
990 escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
991
992 boost::python::list keys=arg.keys();
993 for (int i=0;i<num_data;++i) {
994 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
995 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
996 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
997 data[i]=d.getDataC();
998 ptr_data[i]=&(data[i]);
999 std::string n=boost::python::extract<std::string>(keys[i]);
1000 c_names[i]=&(names[i][0]);
1001 if (n.length()>MAX_namelength-1) {
1002 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
1003 c_names[i][MAX_namelength-1]='\0';
1004 } else {
1005 strcpy(c_names[i],n.c_str());
1006 }
1007 }
1008 #ifndef PASO_MPI
1009 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1010 #else
1011 Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
1012 #endif
1013
1014 checkFinleyError();
1015 /* win32 refactor */
1016 TMPMEMFREE(c_names);
1017 TMPMEMFREE(data);
1018 TMPMEMFREE(ptr_data);
1019 for(int i=0;i<num_data;i++)
1020 {
1021 TMPMEMFREE(names[i]);
1022 }
1023 TMPMEMFREE(names);
1024
1025 return;
1026 }
1027
1028
1029 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1030 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1031 const int row_blocksize,
1032 const escript::FunctionSpace& row_functionspace,
1033 const int column_blocksize,
1034 const escript::FunctionSpace& column_functionspace,
1035 const int type) const
1036 {
1037 int reduceRowOrder=0;
1038 int reduceColOrder=0;
1039 // is the domain right?
1040 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1041 if (row_domain!=*this)
1042 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1043 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1044 if (col_domain!=*this)
1045 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1046 // is the function space type right
1047 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1048 reduceRowOrder=0;
1049 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1050 reduceRowOrder=1;
1051 } else {
1052 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1053 }
1054 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1055 reduceColOrder=0;
1056 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1057 reduceColOrder=1;
1058 } else {
1059 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1060 }
1061 // generate matrix:
1062
1063 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1064 checkFinleyError();
1065 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1066 checkPasoError();
1067 Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
1068 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1069 }
1070
1071 //
1072 // vtkObject MeshAdapter::createVtkObject() const
1073 // TODO:
1074 //
1075
1076 //
1077 // returns true if data at the atom_type is considered as being cell centered:
1078 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1079 {
1080 switch(functionSpaceCode) {
1081 case(Nodes):
1082 case(DegreesOfFreedom):
1083 case(ReducedDegreesOfFreedom):
1084 return false;
1085 break;
1086 case(Elements):
1087 case(FaceElements):
1088 case(Points):
1089 case(ContactElementsZero):
1090 case(ContactElementsOne):
1091 case(ReducedElements):
1092 case(ReducedFaceElements):
1093 case(ReducedContactElementsZero):
1094 case(ReducedContactElementsOne):
1095 return true;
1096 break;
1097 default:
1098 stringstream temp;
1099 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1100 throw FinleyAdapterException(temp.str());
1101 break;
1102 }
1103 checkFinleyError();
1104 return false;
1105 }
1106
1107 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1108 {
1109 switch(functionSpaceType_source) {
1110 case(Nodes):
1111 switch(functionSpaceType_target) {
1112 case(Nodes):
1113 case(ReducedNodes):
1114 case(ReducedDegreesOfFreedom):
1115 case(DegreesOfFreedom):
1116 case(Elements):
1117 case(ReducedElements):
1118 case(FaceElements):
1119 case(ReducedFaceElements):
1120 case(Points):
1121 case(ContactElementsZero):
1122 case(ReducedContactElementsZero):
1123 case(ContactElementsOne):
1124 case(ReducedContactElementsOne):
1125 return true;
1126 default:
1127 stringstream temp;
1128 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1129 throw FinleyAdapterException(temp.str());
1130 }
1131 break;
1132 case(ReducedNodes):
1133 switch(functionSpaceType_target) {
1134 case(ReducedNodes):
1135 case(ReducedDegreesOfFreedom):
1136 case(Elements):
1137 case(ReducedElements):
1138 case(FaceElements):
1139 case(ReducedFaceElements):
1140 case(Points):
1141 case(ContactElementsZero):
1142 case(ReducedContactElementsZero):
1143 case(ContactElementsOne):
1144 case(ReducedContactElementsOne):
1145 return true;
1146 case(Nodes):
1147 case(DegreesOfFreedom):
1148 return false;
1149 default:
1150 stringstream temp;
1151 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1152 throw FinleyAdapterException(temp.str());
1153 }
1154 break;
1155 case(Elements):
1156 if (functionSpaceType_target==Elements) {
1157 return true;
1158 } else if (functionSpaceType_target==ReducedElements) {
1159 return true;
1160 } else {
1161 return false;
1162 }
1163 case(ReducedElements):
1164 if (functionSpaceType_target==ReducedElements) {
1165 return true;
1166 } else {
1167 return false;
1168 }
1169 case(FaceElements):
1170 if (functionSpaceType_target==FaceElements) {
1171 return true;
1172 } else if (functionSpaceType_target==ReducedFaceElements) {
1173 return true;
1174 } else {
1175 return false;
1176 }
1177 case(ReducedFaceElements):
1178 if (functionSpaceType_target==ReducedFaceElements) {
1179 return true;
1180 } else {
1181 return false;
1182 }
1183 case(Points):
1184 if (functionSpaceType_target==Points) {
1185 return true;
1186 } else {
1187 return false;
1188 }
1189 case(ContactElementsZero):
1190 case(ContactElementsOne):
1191 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1192 return true;
1193 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1194 return true;
1195 } else {
1196 return false;
1197 }
1198 case(ReducedContactElementsZero):
1199 case(ReducedContactElementsOne):
1200 if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1201 return true;
1202 } else {
1203 return false;
1204 }
1205 case(DegreesOfFreedom):
1206 switch(functionSpaceType_target) {
1207 case(ReducedDegreesOfFreedom):
1208 case(DegreesOfFreedom):
1209 case(Nodes):
1210 case(ReducedNodes):
1211 case(Elements):
1212 case(ReducedElements):
1213 case(Points):
1214 case(FaceElements):
1215 case(ReducedFaceElements):
1216 case(ContactElementsZero):
1217 case(ReducedContactElementsZero):
1218 case(ContactElementsOne):
1219 case(ReducedContactElementsOne):
1220 return true;
1221 default:
1222 stringstream temp;
1223 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1224 throw FinleyAdapterException(temp.str());
1225 }
1226 break;
1227 case(ReducedDegreesOfFreedom):
1228 switch(functionSpaceType_target) {
1229 case(ReducedDegreesOfFreedom):
1230 case(ReducedNodes):
1231 case(Elements):
1232 case(ReducedElements):
1233 case(FaceElements):
1234 case(ReducedFaceElements):
1235 case(Points):
1236 case(ContactElementsZero):
1237 case(ReducedContactElementsZero):
1238 case(ContactElementsOne):
1239 case(ReducedContactElementsOne):
1240 return true;
1241 case(Nodes):
1242 case(DegreesOfFreedom):
1243 return false;
1244 default:
1245 stringstream temp;
1246 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1247 throw FinleyAdapterException(temp.str());
1248 }
1249 break;
1250 default:
1251 stringstream temp;
1252 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1253 throw FinleyAdapterException(temp.str());
1254 break;
1255 }
1256 checkFinleyError();
1257 return false;
1258 }
1259
1260 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1261 {
1262 return false;
1263 }
1264
1265 bool MeshAdapter::operator==(const AbstractDomain& other) const
1266 {
1267 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1268 if (temp!=0) {
1269 return (m_finleyMesh==temp->m_finleyMesh);
1270 } else {
1271 return false;
1272 }
1273 }
1274
1275 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1276 {
1277 return !(operator==(other));
1278 }
1279
1280 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1281 {
1282 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1283 checkPasoError();
1284 return out;
1285 }
1286
1287 escript::Data MeshAdapter::getX() const
1288 {
1289 return continuousFunction(asAbstractContinuousDomain()).getX();
1290 }
1291
1292 escript::Data MeshAdapter::getNormal() const
1293 {
1294 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1295 }
1296
1297 escript::Data MeshAdapter::getSize() const
1298 {
1299 return function(asAbstractContinuousDomain()).getSize();
1300 }
1301
1302 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1303 {
1304 int *out=0,i;
1305 Finley_Mesh* mesh=m_finleyMesh.get();
1306 switch (functionSpaceType) {
1307 case(Nodes):
1308 if (mesh->Nodes!=NULL) {
1309 out=mesh->Nodes->Id;
1310 break;
1311 }
1312 case(ReducedNodes):
1313 throw FinleyAdapterException("Error - ReducedNodes not supported yet.");
1314 break;
1315 case(Elements):
1316 out=mesh->Elements->Id;
1317 break;
1318 case(ReducedElements):
1319 out=mesh->Elements->Id;
1320 break;
1321 case(FaceElements):
1322 out=mesh->FaceElements->Id;
1323 break;
1324 case(ReducedFaceElements):
1325 out=mesh->FaceElements->Id;
1326 break;
1327 case(Points):
1328 out=mesh->Points->Id;
1329 break;
1330 case(ContactElementsZero):
1331 out=mesh->ContactElements->Id;
1332 break;
1333 case(ReducedContactElementsZero):
1334 out=mesh->ContactElements->Id;
1335 break;
1336 case(ContactElementsOne):
1337 out=mesh->ContactElements->Id;
1338 break;
1339 case(ReducedContactElementsOne):
1340 out=mesh->ContactElements->Id;
1341 break;
1342 case(DegreesOfFreedom):
1343 out=mesh->Nodes->degreeOfFreedomId;
1344 break;
1345 case(ReducedDegreesOfFreedom):
1346 out=mesh->Nodes->reducedDegreeOfFreedomId;
1347 break;
1348 default:
1349 stringstream temp;
1350 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1351 throw FinleyAdapterException(temp.str());
1352 break;
1353 }
1354 return out;
1355 }
1356 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1357 {
1358 int out=0;
1359 Finley_Mesh* mesh=m_finleyMesh.get();
1360 switch (functionSpaceType) {
1361 case(Nodes):
1362 out=mesh->Nodes->Tag[sampleNo];
1363 break;
1364 case(ReducedNodes):
1365 throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1366 break;
1367 case(Elements):
1368 out=mesh->Elements->Tag[sampleNo];
1369 break;
1370 case(ReducedElements):
1371 out=mesh->Elements->Tag[sampleNo];
1372 break;
1373 case(FaceElements):
1374 out=mesh->FaceElements->Tag[sampleNo];
1375 break;
1376 case(ReducedFaceElements):
1377 out=mesh->FaceElements->Tag[sampleNo];
1378 break;
1379 case(Points):
1380 out=mesh->Points->Tag[sampleNo];
1381 break;
1382 case(ContactElementsZero):
1383 out=mesh->ContactElements->Tag[sampleNo];
1384 break;
1385 case(ReducedContactElementsZero):
1386 out=mesh->ContactElements->Tag[sampleNo];
1387 break;
1388 case(ContactElementsOne):
1389 out=mesh->ContactElements->Tag[sampleNo];
1390 break;
1391 case(ReducedContactElementsOne):
1392 out=mesh->ContactElements->Tag[sampleNo];
1393 break;
1394 case(DegreesOfFreedom):
1395 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1396 break;
1397 case(ReducedDegreesOfFreedom):
1398 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1399 break;
1400 default:
1401 stringstream temp;
1402 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1403 throw FinleyAdapterException(temp.str());
1404 break;
1405 }
1406 return out;
1407 }
1408
1409
1410 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1411 {
1412 Finley_Mesh* mesh=m_finleyMesh.get();
1413 escriptDataC tmp=mask.getDataC();
1414 switch(functionSpaceType) {
1415 case(Nodes):
1416 Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1417 break;
1418 case(ReducedNodes):
1419 throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1420 break;
1421 case(DegreesOfFreedom):
1422 throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1423 break;
1424 case(ReducedDegreesOfFreedom):
1425 throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1426 break;
1427 case(Elements):
1428 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1429 break;
1430 case(ReducedElements):
1431 Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1432 break;
1433 case(FaceElements):
1434 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1435 break;
1436 case(ReducedFaceElements):
1437 Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1438 break;
1439 case(Points):
1440 Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1441 break;
1442 case(ContactElementsZero):
1443 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1444 break;
1445 case(ReducedContactElementsZero):
1446 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1447 break;
1448 case(ContactElementsOne):
1449 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1450 break;
1451 case(ReducedContactElementsOne):
1452 Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1453 break;
1454 default:
1455 stringstream temp;
1456 temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1457 throw FinleyAdapterException(temp.str());
1458 }
1459 checkFinleyError();
1460 return;
1461 }
1462
1463 void MeshAdapter::setTagMap(const std::string& name, int tag)
1464 {
1465 Finley_Mesh* mesh=m_finleyMesh.get();
1466 Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1467 checkPasoError();
1468 // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1469 }
1470
1471 int MeshAdapter::getTag(const std::string& name) const
1472 {
1473 Finley_Mesh* mesh=m_finleyMesh.get();
1474 int tag=0;
1475 tag=Finley_Mesh_getTag(mesh, name.c_str());
1476 checkPasoError();
1477 // throwStandardException("MeshAdapter::getTag is not implemented.");
1478 return tag;
1479 }
1480
1481 bool MeshAdapter::isValidTagName(const std::string& name) const
1482 {
1483 Finley_Mesh* mesh=m_finleyMesh.get();
1484 return Finley_Mesh_isValidTagName(mesh,name.c_str());
1485 }
1486
1487 std::string MeshAdapter::showTagNames() const
1488 {
1489 stringstream temp;
1490 Finley_Mesh* mesh=m_finleyMesh.get();
1491 Finley_TagMap* tag_map=mesh->TagMap;
1492 while (tag_map) {
1493 temp << tag_map->name;
1494 tag_map=tag_map->next;
1495 if (tag_map) temp << ", ";
1496 }
1497 return string(temp.str());
1498 }
1499
1500 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26