/[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 798 - (show annotations)
Fri Aug 4 01:05:36 2006 UTC (12 years, 10 months ago) by gross
File size: 42180 byte(s)
Reimplementation of the assemblage with persistent jacobeans.
There are also a few changes to the tests which has now
dramatically reduced the memory demand.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26