/[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 967 - (show annotations)
Tue Feb 13 09:40:12 2007 UTC (12 years, 2 months ago) by gross
File size: 41890 byte(s)
dump and load of expanded data via netCDF added. some test are still missing.
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