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

Contents of /trunk-mpi-branch/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26