/[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 751 - (show annotations)
Mon Jun 26 01:46:34 2006 UTC (12 years, 10 months ago) by bcumming
File size: 39985 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26