/[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 971 - (show annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 8 months ago) by ksteube
File size: 41890 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26