/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Contents of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 757 - (show annotations)
Mon Jun 26 13:12:56 2006 UTC (12 years, 10 months ago) by woo409
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 41135 byte(s)
+ Merge of intelc_win32 branch (revision 741:755) with trunk. Tested on iVEC altix (run_tests and py_tests all pass)

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26