/[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 682 - (show annotations)
Mon Mar 27 02:43:09 2006 UTC (13 years ago) by robwdcock
File size: 37983 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


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 numSamples=mesh->Nodes->numDegreesOfFreedom;
220 }
221 break;
222 case(ReducedDegreesOfFreedom):
223 if (mesh->Nodes!=NULL) {
224 numDataPointsPerSample=1;
225 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
226 }
227 break;
228 default:
229 stringstream temp;
230 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
231 throw FinleyAdapterException(temp.str());
232 break;
233 }
234 return pair<int,int>(numDataPointsPerSample,numSamples);
235 }
236
237 //
238 // adds linear PDE of second order into a given stiffness matrix and right hand side:
239 //
240 void MeshAdapter::addPDEToSystem(
241 SystemMatrixAdapter& mat, escript::Data& rhs,
242 const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
243 const escript::Data& d, const escript::Data& y,
244 const escript::Data& d_contact,const escript::Data& y_contact) const
245 {
246 Finley_Mesh* mesh=m_finleyMesh.get();
247 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),
248 &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));
249 checkFinleyError();
250 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,
251 mat.getPaso_SystemMatrix(),
252 &(rhs.getDataC()),
253 &(d.getDataC()),&(y.getDataC()),
254 Finley_Assemble_handelShapeMissMatch_Mean_out);
255 checkFinleyError();
256 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,
257 mat.getPaso_SystemMatrix(),
258 &(rhs.getDataC()),
259 &(d_contact.getDataC()),
260 &(y_contact.getDataC()),
261 Finley_Assemble_handelShapeMissMatch_Step_out);
262 checkFinleyError();
263 }
264
265 //
266 // adds linear PDE of second order into the right hand side only
267 //
268 void MeshAdapter::addPDEToRHS( escript::Data& rhs,
269 const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
270 {
271 Finley_Mesh* mesh=m_finleyMesh.get();
272
273 // Finley_Assemble_PDE_RHS(mesh->Nodes,mesh->Elements,&(rhs.getDataC()),&(X.getDataC()),&(Y.getDataC()));
274 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,0,&(rhs.getDataC()),0,0,0,0,&(X.getDataC()),&(Y.getDataC()));
275 checkFinleyError();
276
277 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->FaceElements,&(rhs.getDataC()),&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
278 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,0,&(rhs.getDataC()),0,&(y.getDataC()),Finley_Assemble_handelShapeMissMatch_Mean_out);
279
280 checkFinleyError();
281 Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,0,&(rhs.getDataC()),0,&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
282 // Finley_Assemble_RobinCondition_RHS(mesh->Nodes,mesh->ContactElements,&(rhs.getDataC()),&(y_contact.getDataC()),Finley_Assemble_handelShapeMissMatch_Step_out);
283 checkFinleyError();
284 }
285
286 //
287 // interpolates data between different function spaces:
288 //
289 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
290 {
291 const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
292 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
293 if (inDomain!=*this)
294 throw FinleyAdapterException("Error - Illegal domain of interpolant.");
295 if (targetDomain!=*this)
296 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
297
298 Finley_Mesh* mesh=m_finleyMesh.get();
299 switch(in.getFunctionSpace().getTypeCode()) {
300 case(Nodes):
301 switch(target.getFunctionSpace().getTypeCode()) {
302 case(Nodes):
303 case(ReducedDegreesOfFreedom):
304 case(DegreesOfFreedom):
305 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
306 break;
307 case(Elements):
308 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
309 break;
310 case(FaceElements):
311 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
312 break;
313 case(Points):
314 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
315 break;
316 case(ContactElementsZero):
317 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
318 break;
319 case(ContactElementsOne):
320 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
321 break;
322 default:
323 stringstream temp;
324 temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
325 throw FinleyAdapterException(temp.str());
326 break;
327 }
328 break;
329 case(Elements):
330 if (target.getFunctionSpace().getTypeCode()==Elements) {
331 Finley_Assemble_CopyElementData(mesh->Elements,&(target.getDataC()),&(in.getDataC()));
332 } else {
333 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
334 }
335 break;
336 case(FaceElements):
337 if (target.getFunctionSpace().getTypeCode()==FaceElements) {
338 Finley_Assemble_CopyElementData(mesh->FaceElements,&(target.getDataC()),&(in.getDataC()));
339 } else {
340 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
341 break;
342 }
343 case(Points):
344 if (target.getFunctionSpace().getTypeCode()==Points) {
345 Finley_Assemble_CopyElementData(mesh->Points,&(target.getDataC()),&(in.getDataC()));
346 } else {
347 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
348 break;
349 }
350 break;
351 case(ContactElementsZero):
352 case(ContactElementsOne):
353 if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
354 Finley_Assemble_CopyElementData(mesh->ContactElements,&(target.getDataC()),&(in.getDataC()));
355 } else {
356 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
357 break;
358 }
359 break;
360 case(DegreesOfFreedom):
361 switch(target.getFunctionSpace().getTypeCode()) {
362 case(ReducedDegreesOfFreedom):
363 case(DegreesOfFreedom):
364 case(Nodes):
365 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
366 break;
367 case(Elements):
368 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
369 break;
370 case(FaceElements):
371 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
372 break;
373 case(Points):
374 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
375 break;
376 case(ContactElementsZero):
377 case(ContactElementsOne):
378 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
379 break;
380 default:
381 stringstream temp;
382 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
383 throw FinleyAdapterException(temp.str());
384 break;
385 }
386 break;
387 case(ReducedDegreesOfFreedom):
388 switch(target.getFunctionSpace().getTypeCode()) {
389 case(ReducedDegreesOfFreedom):
390 Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
391 break;
392 case(Elements):
393 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
394 break;
395 case(FaceElements):
396 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(in.getDataC()),&(target.getDataC()));
397 break;
398 case(Points):
399 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(in.getDataC()),&(target.getDataC()));
400 break;
401 case(ContactElementsZero):
402 case(ContactElementsOne):
403 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
404 break;
405 case(Nodes):
406 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
407 break;
408 case(DegreesOfFreedom):
409 throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
410 break;
411 default:
412 stringstream temp;
413 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
414 throw FinleyAdapterException(temp.str());
415 break;
416 }
417 break;
418 default:
419 stringstream temp;
420 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
421 throw FinleyAdapterException(temp.str());
422 break;
423 }
424 checkFinleyError();
425 }
426
427 //
428 // copies the locations of sample points into x:
429 //
430 void MeshAdapter::setToX(escript::Data& arg) const
431 {
432 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
433 if (argDomain!=*this)
434 throw FinleyAdapterException("Error - Illegal domain of data point locations");
435 Finley_Mesh* mesh=m_finleyMesh.get();
436 // in case of values node coordinates we can do the job directly:
437 if (arg.getFunctionSpace().getTypeCode()==Nodes) {
438 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(arg.getDataC()));
439 } else {
440 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
441 Finley_Assemble_NodeCoordinates(mesh->Nodes,&(tmp_data.getDataC()));
442 // this is then interpolated onto arg:
443 interpolateOnDomain(arg,tmp_data);
444 }
445 checkFinleyError();
446 }
447
448 //
449 // return the normal vectors at the location of data points as a Data object:
450 //
451 void MeshAdapter::setToNormal(escript::Data& normal) const
452 {
453 const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
454 if (normalDomain!=*this)
455 throw FinleyAdapterException("Error - Illegal domain of normal locations");
456 Finley_Mesh* mesh=m_finleyMesh.get();
457 switch(normal.getFunctionSpace().getTypeCode()) {
458 case(Nodes):
459 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
460 break;
461 case(Elements):
462 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
463 break;
464 case (FaceElements):
465 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&(normal.getDataC()));
466 break;
467 case(Points):
468 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
469 break;
470 case (ContactElementsOne):
471 case (ContactElementsZero):
472 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&(normal.getDataC()));
473 break;
474 case(DegreesOfFreedom):
475 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
476 break;
477 case(ReducedDegreesOfFreedom):
478 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
479 break;
480 default:
481 stringstream temp;
482 temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
483 throw FinleyAdapterException(temp.str());
484 break;
485 }
486 checkFinleyError();
487 }
488
489 //
490 // interpolates data to other domain:
491 //
492 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
493 {
494 const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
495 if (targetDomain!=*this)
496 throw FinleyAdapterException("Error - Illegal domain of interpolation target");
497
498 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
499 }
500
501 //
502 // calculates the integral of a function defined of arg:
503 //
504 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
505 {
506 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
507 if (argDomain!=*this)
508 throw FinleyAdapterException("Error - Illegal domain of integration kernel");
509
510 Finley_Mesh* mesh=m_finleyMesh.get();
511 switch(arg.getFunctionSpace().getTypeCode()) {
512 case(Nodes):
513 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
514 break;
515 case(Elements):
516 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&(arg.getDataC()),&integrals[0]);
517 break;
518 case(FaceElements):
519 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&(arg.getDataC()),&integrals[0]);
520 break;
521 case(Points):
522 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
523 break;
524 case(ContactElementsZero):
525 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
526 break;
527 case(ContactElementsOne):
528 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&(arg.getDataC()),&integrals[0]);
529 break;
530 case(DegreesOfFreedom):
531 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
532 break;
533 case(ReducedDegreesOfFreedom):
534 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
535 break;
536 default:
537 stringstream temp;
538 temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
539 throw FinleyAdapterException(temp.str());
540 break;
541 }
542 checkFinleyError();
543 }
544
545 //
546 // calculates the gradient of arg:
547 //
548 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
549 {
550 const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
551 if (argDomain!=*this)
552 throw FinleyAdapterException("Error - Illegal domain of gradient argument");
553 const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
554 if (gradDomain!=*this)
555 throw FinleyAdapterException("Error - Illegal domain of gradient");
556
557 Finley_Mesh* mesh=m_finleyMesh.get();
558 switch(grad.getFunctionSpace().getTypeCode()) {
559 case(Nodes):
560 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
561 break;
562 case(Elements):
563 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));
564 break;
565 case(FaceElements):
566 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));
567 break;
568 case(Points):
569 throw FinleyAdapterException("Error - Gradient at points is not supported.");
570 break;
571 case(ContactElementsZero):
572 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
573 break;
574 case(ContactElementsOne):
575 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));
576 break;
577 case(DegreesOfFreedom):
578 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
579 break;
580 case(ReducedDegreesOfFreedom):
581 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
582 break;
583 default:
584 stringstream temp;
585 temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
586 throw FinleyAdapterException(temp.str());
587 break;
588 }
589 checkFinleyError();
590 }
591
592 //
593 // returns the size of elements:
594 //
595 void MeshAdapter::setToSize(escript::Data& size) const
596 {
597 Finley_Mesh* mesh=m_finleyMesh.get();
598 escriptDataC tmp=size.getDataC();
599 switch(size.getFunctionSpace().getTypeCode()) {
600 case(Nodes):
601 throw FinleyAdapterException("Error - Size of nodes is not supported.");
602 break;
603 case(Elements):
604 Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
605 break;
606 case(FaceElements):
607 Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
608 break;
609 case(Points):
610 throw FinleyAdapterException("Error - Size of point elements is not supported.");
611 break;
612 case(ContactElementsZero):
613 case(ContactElementsOne):
614 Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
615 break;
616 case(DegreesOfFreedom):
617 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
618 break;
619 case(ReducedDegreesOfFreedom):
620 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
621 break;
622 default:
623 stringstream temp;
624 temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
625 throw FinleyAdapterException(temp.str());
626 break;
627 }
628 checkFinleyError();
629 }
630
631 // sets the location of nodes:
632 void MeshAdapter::setNewX(const escript::Data& new_x)
633 {
634 Finley_Mesh* mesh=m_finleyMesh.get();
635 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
636 if (newDomain!=*this)
637 throw FinleyAdapterException("Error - Illegal domain of new point locations");
638 Finley_Mesh_setCoordinates(mesh,&(new_x.getDataC()));
639 checkFinleyError();
640 }
641
642 // saves a data array in openDX format:
643 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
644 {
645 int MAX_namelength=256;
646 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
647 char names[num_data][MAX_namelength];
648 char* c_names[num_data];
649 escriptDataC data[num_data];
650 escriptDataC* ptr_data[num_data];
651
652 boost::python::list keys=arg.keys();
653 for (int i=0;i<num_data;++i) {
654 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
655 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
656 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
657 data[i]=d.getDataC();
658 ptr_data[i]=&(data[i]);
659 std::string n=boost::python::extract<std::string>(keys[i]);
660 c_names[i]=&(names[i][0]);
661 if (n.length()>MAX_namelength-1) {
662 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
663 c_names[i][MAX_namelength-1]='\0';
664 } else {
665 strcpy(c_names[i],n.c_str());
666 }
667 }
668 Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
669 checkFinleyError();
670 return;
671 }
672
673 // saves a data array in openVTK format:
674 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
675 {
676 int MAX_namelength=256;
677 const int num_data=boost::python::extract<int>(arg.attr("__len__")());
678 char names[num_data][MAX_namelength];
679 char* c_names[num_data];
680 escriptDataC data[num_data];
681 escriptDataC* ptr_data[num_data];
682
683 boost::python::list keys=arg.keys();
684 for (int i=0;i<num_data;++i) {
685 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
686 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
687 throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
688 data[i]=d.getDataC();
689 ptr_data[i]=&(data[i]);
690 std::string n=boost::python::extract<std::string>(keys[i]);
691 c_names[i]=&(names[i][0]);
692 if (n.length()>MAX_namelength-1) {
693 strncpy(c_names[i],n.c_str(),MAX_namelength-1);
694 c_names[i][MAX_namelength-1]='\0';
695 } else {
696 strcpy(c_names[i],n.c_str());
697 }
698 }
699 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
700 checkFinleyError();
701 return;
702 }
703
704
705 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
706 SystemMatrixAdapter MeshAdapter::newSystemMatrix(
707 const int row_blocksize,
708 const escript::FunctionSpace& row_functionspace,
709 const int column_blocksize,
710 const escript::FunctionSpace& column_functionspace,
711 const int type) const
712 {
713 int reduceRowOrder=0;
714 int reduceColOrder=0;
715 // is the domain right?
716 const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
717 if (row_domain!=*this)
718 throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
719 const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
720 if (col_domain!=*this)
721 throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
722 // is the function space type right
723 if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
724 reduceRowOrder=0;
725 } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
726 reduceRowOrder=1;
727 } else {
728 throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
729 }
730 if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
731 reduceColOrder=0;
732 } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
733 reduceColOrder=1;
734 } else {
735 throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
736 }
737 // generate matrix:
738
739 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
740 checkFinleyError();
741 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
742 checkPasoError();
743 Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
744 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
745 }
746
747 //
748 // vtkObject MeshAdapter::createVtkObject() const
749 // TODO:
750 //
751
752 //
753 // returns true if data at the atom_type is considered as being cell centered:
754 bool MeshAdapter::isCellOriented(int functionSpaceCode) const
755 {
756 switch(functionSpaceCode) {
757 case(Nodes):
758 case(DegreesOfFreedom):
759 case(ReducedDegreesOfFreedom):
760 return false;
761 break;
762 case(Elements):
763 case(FaceElements):
764 case(Points):
765 case(ContactElementsZero):
766 case(ContactElementsOne):
767 return true;
768 break;
769 default:
770 stringstream temp;
771 temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
772 throw FinleyAdapterException(temp.str());
773 break;
774 }
775 checkFinleyError();
776 return false;
777 }
778
779 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
780 {
781 switch(functionSpaceType_source) {
782 case(Nodes):
783 switch(functionSpaceType_target) {
784 case(Nodes):
785 case(ReducedDegreesOfFreedom):
786 case(DegreesOfFreedom):
787 case(Elements):
788 case(FaceElements):
789 case(Points):
790 case(ContactElementsZero):
791 case(ContactElementsOne):
792 return true;
793 default:
794 stringstream temp;
795 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
796 throw FinleyAdapterException(temp.str());
797 }
798 break;
799 case(Elements):
800 if (functionSpaceType_target==Elements) {
801 return true;
802 } else {
803 return false;
804 }
805 case(FaceElements):
806 if (functionSpaceType_target==FaceElements) {
807 return true;
808 } else {
809 return false;
810 }
811 case(Points):
812 if (functionSpaceType_target==Points) {
813 return true;
814 } else {
815 return false;
816 }
817 case(ContactElementsZero):
818 case(ContactElementsOne):
819 if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
820 return true;
821 } else {
822 return false;
823 }
824 case(DegreesOfFreedom):
825 switch(functionSpaceType_target) {
826 case(ReducedDegreesOfFreedom):
827 case(DegreesOfFreedom):
828 case(Nodes):
829 case(Elements):
830 case(FaceElements):
831 case(Points):
832 case(ContactElementsZero):
833 case(ContactElementsOne):
834 return true;
835 default:
836 stringstream temp;
837 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
838 throw FinleyAdapterException(temp.str());
839 }
840 break;
841 case(ReducedDegreesOfFreedom):
842 switch(functionSpaceType_target) {
843 case(ReducedDegreesOfFreedom):
844 case(Elements):
845 case(FaceElements):
846 case(Points):
847 case(ContactElementsZero):
848 case(ContactElementsOne):
849 return true;
850 case(Nodes):
851 case(DegreesOfFreedom):
852 return false;
853 default:
854 stringstream temp;
855 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
856 throw FinleyAdapterException(temp.str());
857 }
858 break;
859 default:
860 stringstream temp;
861 temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
862 throw FinleyAdapterException(temp.str());
863 break;
864 }
865 checkFinleyError();
866 return false;
867 }
868
869 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
870 {
871 return false;
872 }
873
874 bool MeshAdapter::operator==(const AbstractDomain& other) const
875 {
876 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
877 if (temp!=0) {
878 return (m_finleyMesh==temp->m_finleyMesh);
879 } else {
880 return false;
881 }
882 }
883
884 bool MeshAdapter::operator!=(const AbstractDomain& other) const
885 {
886 return !(operator==(other));
887 }
888
889 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
890 {
891 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
892 checkPasoError();
893 return out;
894 }
895
896 escript::Data MeshAdapter::getX() const
897 {
898 return continuousFunction(asAbstractContinuousDomain()).getX();
899 }
900
901 escript::Data MeshAdapter::getNormal() const
902 {
903 return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
904 }
905
906 escript::Data MeshAdapter::getSize() const
907 {
908 return function(asAbstractContinuousDomain()).getSize();
909 }
910
911 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
912 {
913 int out=0;
914 Finley_Mesh* mesh=m_finleyMesh.get();
915 switch (functionSpaceType) {
916 case(Nodes):
917 out=mesh->Nodes->Tag[sampleNo];
918 break;
919 case(Elements):
920 out=mesh->Elements->Tag[sampleNo];
921 break;
922 case(FaceElements):
923 out=mesh->FaceElements->Tag[sampleNo];
924 break;
925 case(Points):
926 out=mesh->Points->Tag[sampleNo];
927 break;
928 case(ContactElementsZero):
929 out=mesh->ContactElements->Tag[sampleNo];
930 break;
931 case(ContactElementsOne):
932 out=mesh->ContactElements->Tag[sampleNo];
933 break;
934 case(DegreesOfFreedom):
935 break;
936 case(ReducedDegreesOfFreedom):
937 break;
938 default:
939 stringstream temp;
940 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
941 throw FinleyAdapterException(temp.str());
942 break;
943 }
944 return out;
945 }
946 int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
947 {
948 int out=0,i;
949 Finley_Mesh* mesh=m_finleyMesh.get();
950 switch (functionSpaceType) {
951 case(Nodes):
952 if (mesh->Nodes!=NULL) {
953 out=mesh->Nodes->Id[sampleNo];
954 break;
955 }
956 case(Elements):
957 out=mesh->Elements->Id[sampleNo];
958 break;
959 case(FaceElements):
960 out=mesh->FaceElements->Id[sampleNo];
961 break;
962 case(Points):
963 out=mesh->Points->Id[sampleNo];
964 break;
965 case(ContactElementsZero):
966 out=mesh->ContactElements->Id[sampleNo];
967 break;
968 case(ContactElementsOne):
969 out=mesh->ContactElements->Id[sampleNo];
970 break;
971 case(DegreesOfFreedom):
972 for (i=0;i<mesh->Nodes->numNodes; ++i) {
973 if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
974 out=mesh->Nodes->Id[i];
975 break;
976 }
977 }
978 break;
979 case(ReducedDegreesOfFreedom):
980 for (i=0;i<mesh->Nodes->numNodes; ++i) {
981 if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
982 out=mesh->Nodes->Id[i];
983 break;
984 }
985 }
986 break;
987 default:
988 stringstream temp;
989 temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
990 throw FinleyAdapterException(temp.str());
991 break;
992 }
993 return out;
994 }
995
996 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26