/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Annotation of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1137 - (hide annotations)
Thu May 10 08:11:31 2007 UTC (11 years, 11 months ago) by gross
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 56150 byte(s)
This version passes the tests on windows except for 

   * vtk
   * netCDF

The version needs to be tested on altix and linux
1 jgs 102 // $Id$
2 jgs 82 /*
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 jgs 472
16 ksteube 817 #ifdef PASO_MPI
17     #include <mpi.h>
18     #endif
19 jgs 203 #include "MeshAdapter.h"
20 jgs 149
21 robwdcock 682 #include "escript/Data.h"
22     #include "escript/DataFactory.h"
23 jgs 480
24 jgs 82 using namespace std;
25     using namespace escript;
26    
27     namespace finley {
28 jgs 149
29 jgs 82 //
30 jgs 149 // define the static constants
31 jgs 82 MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
32     const int MeshAdapter::DegreesOfFreedom=FINLEY_DEGREES_OF_FREEDOM;
33     const int MeshAdapter::ReducedDegreesOfFreedom=FINLEY_REDUCED_DEGREES_OF_FREEDOM;
34     const int MeshAdapter::Nodes=FINLEY_NODES;
35 gross 1062 const int MeshAdapter::ReducedNodes=FINLEY_REDUCED_NODES;
36 jgs 82 const int MeshAdapter::Elements=FINLEY_ELEMENTS;
37 gross 1059 const int MeshAdapter::ReducedElements=FINLEY_REDUCED_ELEMENTS;
38 jgs 82 const int MeshAdapter::FaceElements=FINLEY_FACE_ELEMENTS;
39 gross 1059 const int MeshAdapter::ReducedFaceElements=FINLEY_REDUCED_FACE_ELEMENTS;
40 jgs 82 const int MeshAdapter::Points=FINLEY_POINTS;
41     const int MeshAdapter::ContactElementsZero=FINLEY_CONTACT_ELEMENTS_1;
42 gross 1059 const int MeshAdapter::ReducedContactElementsZero=FINLEY_REDUCED_CONTACT_ELEMENTS_1;
43 jgs 82 const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
44 gross 1059 const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
45 jgs 82
46     MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)
47     {
48     setFunctionSpaceTypeNames();
49     //
50     // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer
51     // for us.
52     m_finleyMesh.reset(finleyMesh,null_deleter());
53     }
54 jgs 149
55 jgs 82 //
56     // The copy constructor should just increment the use count
57     MeshAdapter::MeshAdapter(const MeshAdapter& in):
58     m_finleyMesh(in.m_finleyMesh)
59     {
60     setFunctionSpaceTypeNames();
61     }
62    
63     MeshAdapter::~MeshAdapter()
64     {
65     //
66     // I hope the case for the pointer being zero has been taken care of.
67     // cout << "In MeshAdapter destructor." << endl;
68     if (m_finleyMesh.unique()) {
69     // cout << "Calling dealloc." << endl;
70     Finley_Mesh_dealloc(m_finleyMesh.get());
71     // cout << "Finished dealloc." << endl;
72     }
73     }
74    
75     Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
76     return m_finleyMesh.get();
77     }
78    
79     void MeshAdapter::write(const std::string& fileName) const
80     {
81 woo409 757 char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
82 jgs 82 strcpy(fName,fileName.c_str());
83     Finley_Mesh_write(m_finleyMesh.get(),fName);
84     checkFinleyError();
85 woo409 757 TMPMEMFREE(fName);
86 jgs 82 }
87    
88     string MeshAdapter::getDescription() const
89     {
90 gross 1137 return string("FinleyMesh");
91 jgs 82 }
92    
93     string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
94     {
95     FunctionSpaceNamesMapType::iterator loc;
96     loc=m_functionSpaceTypeNames.find(functionSpaceType);
97     if (loc==m_functionSpaceTypeNames.end()) {
98     return "Invalid function space type code.";
99     } else {
100 gross 1137 return string(loc->second);
101 jgs 82 }
102     }
103    
104     bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
105     {
106     FunctionSpaceNamesMapType::iterator loc;
107     loc=m_functionSpaceTypeNames.find(functionSpaceType);
108     return (loc!=m_functionSpaceTypeNames.end());
109     }
110    
111     void MeshAdapter::setFunctionSpaceTypeNames()
112     {
113     m_functionSpaceTypeNames.insert
114     (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Finley_DegreesOfFreedom"));
115     m_functionSpaceTypeNames.insert
116     (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Finley_ReducedDegreesOfFreedom"));
117     m_functionSpaceTypeNames.insert
118     (FunctionSpaceNamesMapType::value_type(Nodes,"Finley_Nodes"));
119     m_functionSpaceTypeNames.insert
120 gross 1062 (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Finley_Reduced_Nodes"));
121     m_functionSpaceTypeNames.insert
122 jgs 82 (FunctionSpaceNamesMapType::value_type(Elements,"Finley_Elements"));
123     m_functionSpaceTypeNames.insert
124 gross 1059 (FunctionSpaceNamesMapType::value_type(ReducedElements,"Finley_Reduced_Elements"));
125     m_functionSpaceTypeNames.insert
126 jgs 82 (FunctionSpaceNamesMapType::value_type(FaceElements,"Finley_Face_Elements"));
127     m_functionSpaceTypeNames.insert
128 gross 1059 (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Finley_Reduced_Face_Elements"));
129     m_functionSpaceTypeNames.insert
130 jgs 82 (FunctionSpaceNamesMapType::value_type(Points,"Finley_Points"));
131     m_functionSpaceTypeNames.insert
132     (FunctionSpaceNamesMapType::value_type(ContactElementsZero,"Finley_Contact_Elements_0"));
133     m_functionSpaceTypeNames.insert
134 gross 1059 (FunctionSpaceNamesMapType::value_type(ReducedContactElementsZero,"Finley_Reduced_Contact_Elements_0"));
135     m_functionSpaceTypeNames.insert
136 jgs 82 (FunctionSpaceNamesMapType::value_type(ContactElementsOne,"Finley_Contact_Elements_1"));
137 gross 1059 m_functionSpaceTypeNames.insert
138     (FunctionSpaceNamesMapType::value_type(ReducedContactElementsOne,"Finley_Reduced_Contact_Elements_1"));
139 jgs 82 }
140    
141     int MeshAdapter::getContinuousFunctionCode() const
142     {
143     return Nodes;
144     }
145 gross 1062 int MeshAdapter::getReducedContinuousFunctionCode() const
146     {
147     return ReducedNodes;
148     }
149 jgs 149
150 jgs 82 int MeshAdapter::getFunctionCode() const
151     {
152     return Elements;
153     }
154 gross 1059 int MeshAdapter::getReducedFunctionCode() const
155     {
156     return ReducedElements;
157     }
158 jgs 149
159 jgs 82 int MeshAdapter::getFunctionOnBoundaryCode() const
160     {
161     return FaceElements;
162     }
163 gross 1059 int MeshAdapter::getReducedFunctionOnBoundaryCode() const
164     {
165     return ReducedFaceElements;
166     }
167 jgs 149
168 jgs 82 int MeshAdapter::getFunctionOnContactZeroCode() const
169     {
170     return ContactElementsZero;
171     }
172 gross 1059 int MeshAdapter::getReducedFunctionOnContactZeroCode() const
173     {
174     return ReducedContactElementsZero;
175     }
176 jgs 149
177 jgs 82 int MeshAdapter::getFunctionOnContactOneCode() const
178     {
179     return ContactElementsOne;
180     }
181 gross 1059 int MeshAdapter::getReducedFunctionOnContactOneCode() const
182     {
183     return ReducedContactElementsOne;
184     }
185 jgs 82
186     int MeshAdapter::getSolutionCode() const
187     {
188     return DegreesOfFreedom;
189     }
190 jgs 149
191 jgs 82 int MeshAdapter::getReducedSolutionCode() const
192     {
193     return ReducedDegreesOfFreedom;
194     }
195 jgs 149
196 jgs 82 int MeshAdapter::getDiracDeltaFunctionCode() const
197     {
198     return Points;
199     }
200 jgs 149
201 jgs 82 //
202     // return the spatial dimension of the Mesh:
203     //
204     int MeshAdapter::getDim() const
205     {
206     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());
207     checkFinleyError();
208     return numDim;
209     }
210 jgs 149
211 jgs 82 //
212     // return the number of data points per sample and the number of samples
213     // needed to represent data on a parts of the mesh.
214     //
215     pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
216     {
217     int numDataPointsPerSample=0;
218     int numSamples=0;
219     Finley_Mesh* mesh=m_finleyMesh.get();
220     switch (functionSpaceCode) {
221     case(Nodes):
222     numDataPointsPerSample=1;
223     if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
224     break;
225 gross 1062 case(ReducedNodes):
226     /* TODO: add ReducedNodes */
227     numDataPointsPerSample=1;
228     if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;
229     throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");
230     break;
231 jgs 82 case(Elements):
232     if (mesh->Elements!=NULL) {
233     numSamples=mesh->Elements->numElements;
234     numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;
235     }
236     break;
237 gross 1062 case(ReducedElements):
238     if (mesh->Elements!=NULL) {
239     numSamples=mesh->Elements->numElements;
240     numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;
241     }
242     break;
243 jgs 82 case(FaceElements):
244     if (mesh->FaceElements!=NULL) {
245     numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;
246     numSamples=mesh->FaceElements->numElements;
247     }
248     break;
249 gross 1062 case(ReducedFaceElements):
250     if (mesh->FaceElements!=NULL) {
251     numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;
252     numSamples=mesh->FaceElements->numElements;
253     }
254     break;
255 jgs 82 case(Points):
256     if (mesh->Points!=NULL) {
257     numDataPointsPerSample=1;
258     numSamples=mesh->Points->numElements;
259     }
260     break;
261     case(ContactElementsZero):
262     if (mesh->ContactElements!=NULL) {
263     numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
264     numSamples=mesh->ContactElements->numElements;
265     }
266     break;
267 gross 1062 case(ReducedContactElementsZero):
268     if (mesh->ContactElements!=NULL) {
269     numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
270     numSamples=mesh->ContactElements->numElements;
271     }
272     break;
273 jgs 82 case(ContactElementsOne):
274     if (mesh->ContactElements!=NULL) {
275     numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;
276     numSamples=mesh->ContactElements->numElements;
277     }
278     break;
279 gross 1062 case(ReducedContactElementsOne):
280     if (mesh->ContactElements!=NULL) {
281     numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;
282     numSamples=mesh->ContactElements->numElements;
283     }
284     break;
285 jgs 82 case(DegreesOfFreedom):
286     if (mesh->Nodes!=NULL) {
287     numDataPointsPerSample=1;
288 bcumming 751 #ifndef PASO_MPI
289 jgs 82 numSamples=mesh->Nodes->numDegreesOfFreedom;
290 bcumming 751 #else
291     numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
292     #endif
293 jgs 82 }
294     break;
295     case(ReducedDegreesOfFreedom):
296     if (mesh->Nodes!=NULL) {
297     numDataPointsPerSample=1;
298 bcumming 751 #ifndef PASO_MPI
299 jgs 82 numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
300 bcumming 751 #else
301     numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
302     #endif
303 jgs 82 }
304     break;
305     default:
306     stringstream temp;
307 jgs 150 temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
308 jgs 82 throw FinleyAdapterException(temp.str());
309     break;
310     }
311     return pair<int,int>(numDataPointsPerSample,numSamples);
312     }
313 jgs 149
314 jgs 82 //
315     // adds linear PDE of second order into a given stiffness matrix and right hand side:
316     //
317     void MeshAdapter::addPDEToSystem(
318 jgs 480 SystemMatrixAdapter& mat, escript::Data& rhs,
319     const escript::Data& A, const escript::Data& B, const escript::Data& C,const escript::Data& D,const escript::Data& X,const escript::Data& Y,
320     const escript::Data& d, const escript::Data& y,
321     const escript::Data& d_contact,const escript::Data& y_contact) const
322 jgs 82 {
323 gross 798 escriptDataC _rhs=rhs.getDataC();
324     escriptDataC _A =A.getDataC();
325     escriptDataC _B=B.getDataC();
326     escriptDataC _C=C.getDataC();
327     escriptDataC _D=D.getDataC();
328     escriptDataC _X=X.getDataC();
329     escriptDataC _Y=Y.getDataC();
330     escriptDataC _d=d.getDataC();
331     escriptDataC _y=y.getDataC();
332     escriptDataC _d_contact=d_contact.getDataC();
333     escriptDataC _y_contact=y_contact.getDataC();
334 bcumming 751
335 jgs 82 Finley_Mesh* mesh=m_finleyMesh.get();
336 gross 798
337 bcumming 751 Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
338 jgs 82 checkFinleyError();
339 bcumming 751
340 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
341 jgs 82 checkFinleyError();
342 bcumming 751
343 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
344 jgs 82 checkFinleyError();
345     }
346 jgs 149
347 jgs 82 //
348 jgs 102 // adds linear PDE of second order into the right hand side only
349     //
350 bcumming 751 void MeshAdapter::addPDEToRHS( escript::Data& rhs, const escript::Data& X,const escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
351 jgs 102 {
352     Finley_Mesh* mesh=m_finleyMesh.get();
353 jgs 147
354 gross 798 escriptDataC _rhs=rhs.getDataC();
355     escriptDataC _X=X.getDataC();
356     escriptDataC _Y=Y.getDataC();
357     escriptDataC _y=y.getDataC();
358     escriptDataC _y_contact=y_contact.getDataC();
359    
360     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
361 jgs 102 checkFinleyError();
362 jgs 147
363 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
364     checkFinleyError();
365 jgs 147
366 gross 798 Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
367 jgs 102 checkFinleyError();
368     }
369 jgs 149
370 jgs 102 //
371 jgs 82 // interpolates data between different function spaces:
372     //
373 jgs 480 void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
374 jgs 82 {
375     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
376     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
377 bcumming 751 if (inDomain!=*this)
378     throw FinleyAdapterException("Error - Illegal domain of interpolant.");
379 jgs 82 if (targetDomain!=*this)
380 bcumming 751 throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
381 jgs 82
382     Finley_Mesh* mesh=m_finleyMesh.get();
383 gross 1062 escriptDataC _target=target.getDataC();
384     escriptDataC _in=in.getDataC();
385 jgs 82 switch(in.getFunctionSpace().getTypeCode()) {
386     case(Nodes):
387     switch(target.getFunctionSpace().getTypeCode()) {
388     case(Nodes):
389 gross 1062 case(ReducedNodes):
390     case(DegreesOfFreedom):
391 jgs 82 case(ReducedDegreesOfFreedom):
392 gross 1062 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
393     break;
394     case(Elements):
395     case(ReducedElements):
396     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
397     break;
398     case(FaceElements):
399     case(ReducedFaceElements):
400     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
401     break;
402     case(Points):
403     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
404     break;
405     case(ContactElementsZero):
406     case(ReducedContactElementsZero):
407     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
408     break;
409     case(ContactElementsOne):
410     case(ReducedContactElementsOne):
411     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
412     break;
413     default:
414     stringstream temp;
415     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
416     throw FinleyAdapterException(temp.str());
417     break;
418     }
419     break;
420     case(ReducedNodes):
421     switch(target.getFunctionSpace().getTypeCode()) {
422     case(Nodes):
423     case(ReducedNodes):
424 jgs 82 case(DegreesOfFreedom):
425 gross 1062 case(ReducedDegreesOfFreedom):
426     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
427 jgs 82 break;
428     case(Elements):
429 gross 1062 case(ReducedElements):
430     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
431 jgs 82 break;
432     case(FaceElements):
433 gross 1062 case(ReducedFaceElements):
434     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
435 jgs 82 break;
436     case(Points):
437 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
438 jgs 82 break;
439     case(ContactElementsZero):
440 gross 1062 case(ReducedContactElementsZero):
441     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
442 jgs 82 break;
443     case(ContactElementsOne):
444 gross 1062 case(ReducedContactElementsOne):
445     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
446 jgs 82 break;
447     default:
448 jgs 150 stringstream temp;
449     temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
450     throw FinleyAdapterException(temp.str());
451 jgs 82 break;
452     }
453     break;
454     case(Elements):
455     if (target.getFunctionSpace().getTypeCode()==Elements) {
456 gross 1062 Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
457 gross 1116 } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
458     Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
459 jgs 82 } else {
460 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
461 jgs 82 }
462     break;
463 gross 1062 case(ReducedElements):
464     if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
465     Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
466     } else {
467     throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
468     }
469     break;
470 jgs 82 case(FaceElements):
471     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
472 gross 1062 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
473 gross 1116 } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
474     Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
475 jgs 82 } else {
476 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
477 jgs 82 }
478 gross 1062 break;
479     case(ReducedFaceElements):
480     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
481     Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
482     } else {
483     throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
484     }
485     break;
486 jgs 82 case(Points):
487     if (target.getFunctionSpace().getTypeCode()==Points) {
488 gross 1062 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
489 jgs 82 } else {
490 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
491 jgs 82 }
492     break;
493     case(ContactElementsZero):
494     case(ContactElementsOne):
495     if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
496 gross 1062 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
497 gross 1116 } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
498     Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
499 jgs 82 } else {
500 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
501 jgs 82 break;
502     }
503     break;
504 gross 1062 case(ReducedContactElementsZero):
505     case(ReducedContactElementsOne):
506     if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
507     Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
508     } else {
509     throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
510     break;
511     }
512     break;
513 bcumming 751 case(DegreesOfFreedom):
514 jgs 82 switch(target.getFunctionSpace().getTypeCode()) {
515     case(ReducedDegreesOfFreedom):
516     case(DegreesOfFreedom):
517     case(Nodes):
518 gross 1062 case(ReducedNodes):
519     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
520 jgs 82 break;
521 bcumming 751 #ifndef PASO_MPI
522 jgs 82 case(Elements):
523 gross 1062 case(ReducedElements):
524     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
525 jgs 82 break;
526     case(FaceElements):
527 gross 1062 case(ReducedFaceElements):
528     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
529 jgs 82 break;
530     case(Points):
531 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
532 jgs 82 break;
533     case(ContactElementsZero):
534     case(ContactElementsOne):
535 gross 1062 case(ReducedContactElementsZero):
536     case(ReducedContactElementsOne):
537     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
538 jgs 82 break;
539 bcumming 751 #else
540     /* need to copy Degrees of freedom data to nodal data so that the external values are available */
541     case(Elements):
542 gross 1062 case(ReducedElements):
543 bcumming 751 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
544 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);
545 bcumming 751 break;
546     case(FaceElements):
547 gross 1062 case(ReducedFaceElements):
548 bcumming 751 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
549 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);
550 bcumming 751 break;
551     case(Points):
552     escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
553 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);
554 bcumming 751 break;
555     case(ContactElementsZero):
556     case(ContactElementsOne):
557 gross 1062 case(ReducedContactElementsZero):
558     case(ReducedContactElementsOne):
559 bcumming 751 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
560 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);
561 bcumming 751 break;
562     #endif
563 jgs 82 default:
564 jgs 150 stringstream temp;
565     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
566     throw FinleyAdapterException(temp.str());
567 jgs 82 break;
568     }
569     break;
570     case(ReducedDegreesOfFreedom):
571     switch(target.getFunctionSpace().getTypeCode()) {
572 gross 1062 case(Nodes):
573     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
574     break;
575     case(ReducedNodes):
576     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
577     break;
578     case(DegreesOfFreedom):
579     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
580     break;
581 jgs 82 case(ReducedDegreesOfFreedom):
582 gross 1062 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
583 jgs 82 break;
584     case(Elements):
585 gross 1062 case(ReducedElements):
586     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
587 jgs 82 break;
588     case(FaceElements):
589 gross 1062 case(ReducedFaceElements):
590     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
591 jgs 82 break;
592     case(Points):
593 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
594 jgs 82 break;
595     case(ContactElementsZero):
596     case(ContactElementsOne):
597 gross 1062 case(ReducedContactElementsZero):
598     case(ReducedContactElementsOne):
599     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
600 jgs 82 break;
601     default:
602 jgs 150 stringstream temp;
603     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
604     throw FinleyAdapterException(temp.str());
605 jgs 82 break;
606     }
607     break;
608     default:
609 jgs 150 stringstream temp;
610     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
611     throw FinleyAdapterException(temp.str());
612 jgs 82 break;
613     }
614     checkFinleyError();
615     }
616    
617     //
618     // copies the locations of sample points into x:
619     //
620 jgs 480 void MeshAdapter::setToX(escript::Data& arg) const
621 jgs 82 {
622     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
623     if (argDomain!=*this)
624     throw FinleyAdapterException("Error - Illegal domain of data point locations");
625     Finley_Mesh* mesh=m_finleyMesh.get();
626     // in case of values node coordinates we can do the job directly:
627     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
628 gross 1062 escriptDataC _arg=arg.getDataC();
629     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
630 jgs 82 } else {
631 jgs 480 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
632 gross 1062 escriptDataC _tmp_data=tmp_data.getDataC();
633     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
634 jgs 82 // this is then interpolated onto arg:
635     interpolateOnDomain(arg,tmp_data);
636     }
637     checkFinleyError();
638     }
639 jgs 149
640 jgs 82 //
641     // return the normal vectors at the location of data points as a Data object:
642     //
643 jgs 480 void MeshAdapter::setToNormal(escript::Data& normal) const
644 jgs 82 {
645     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
646     if (normalDomain!=*this)
647     throw FinleyAdapterException("Error - Illegal domain of normal locations");
648     Finley_Mesh* mesh=m_finleyMesh.get();
649 gross 1062 escriptDataC _normal=normal.getDataC();
650 jgs 82 switch(normal.getFunctionSpace().getTypeCode()) {
651     case(Nodes):
652 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
653 jgs 82 break;
654 gross 1062 case(ReducedNodes):
655     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
656     break;
657 jgs 82 case(Elements):
658 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
659 jgs 82 break;
660 gross 1062 case(ReducedElements):
661     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
662     break;
663 jgs 82 case (FaceElements):
664 gross 1062 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
665 jgs 82 break;
666 gross 1062 case (ReducedFaceElements):
667     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
668     break;
669 jgs 82 case(Points):
670 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
671 jgs 82 break;
672     case (ContactElementsOne):
673     case (ContactElementsZero):
674 gross 1062 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
675 jgs 82 break;
676 gross 1062 case (ReducedContactElementsOne):
677     case (ReducedContactElementsZero):
678     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
679     break;
680 jgs 82 case(DegreesOfFreedom):
681 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
682 jgs 82 break;
683     case(ReducedDegreesOfFreedom):
684 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
685 jgs 82 break;
686     default:
687 jgs 150 stringstream temp;
688     temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
689     throw FinleyAdapterException(temp.str());
690 jgs 82 break;
691     }
692     checkFinleyError();
693     }
694 jgs 149
695 jgs 82 //
696     // interpolates data to other domain:
697     //
698 jgs 480 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
699 jgs 82 {
700     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
701     if (targetDomain!=*this)
702     throw FinleyAdapterException("Error - Illegal domain of interpolation target");
703    
704 jgs 150 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
705 jgs 82 }
706 jgs 149
707 jgs 82 //
708     // calculates the integral of a function defined of arg:
709     //
710 jgs 480 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
711 jgs 82 {
712     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
713     if (argDomain!=*this)
714     throw FinleyAdapterException("Error - Illegal domain of integration kernel");
715    
716     Finley_Mesh* mesh=m_finleyMesh.get();
717 gross 1062 escriptDataC _arg=arg.getDataC();
718 jgs 82 switch(arg.getFunctionSpace().getTypeCode()) {
719     case(Nodes):
720 gross 1062 /* TODO */
721 jgs 150 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
722 jgs 82 break;
723 gross 1062 case(ReducedNodes):
724     /* TODO */
725     throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");
726     break;
727 jgs 82 case(Elements):
728 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
729 jgs 82 break;
730 gross 1062 case(ReducedElements):
731     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
732 gross 1064 break;
733 jgs 82 case(FaceElements):
734 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
735 jgs 82 break;
736 gross 1062 case(ReducedFaceElements):
737     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
738     break;
739 jgs 82 case(Points):
740 jgs 150 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
741 jgs 82 break;
742     case(ContactElementsZero):
743 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
744 jgs 82 break;
745 gross 1062 case(ReducedContactElementsZero):
746     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
747     break;
748 jgs 82 case(ContactElementsOne):
749 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
750 jgs 82 break;
751 gross 1062 case(ReducedContactElementsOne):
752     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
753     break;
754 jgs 82 case(DegreesOfFreedom):
755 jgs 150 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
756 jgs 82 break;
757     case(ReducedDegreesOfFreedom):
758 jgs 150 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
759 jgs 82 break;
760     default:
761 jgs 150 stringstream temp;
762     temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
763     throw FinleyAdapterException(temp.str());
764 jgs 82 break;
765     }
766     checkFinleyError();
767     }
768 jgs 149
769 jgs 82 //
770     // calculates the gradient of arg:
771     //
772 jgs 480 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
773 jgs 82 {
774     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
775     if (argDomain!=*this)
776     throw FinleyAdapterException("Error - Illegal domain of gradient argument");
777     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
778     if (gradDomain!=*this)
779     throw FinleyAdapterException("Error - Illegal domain of gradient");
780    
781     Finley_Mesh* mesh=m_finleyMesh.get();
782 gross 1062 escriptDataC _grad=grad.getDataC();
783 bcumming 751 escriptDataC nodeDataC;
784     #ifdef PASO_MPI
785     escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
786     if( arg.getFunctionSpace().getTypeCode() != Nodes )
787     {
788     Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
789     nodeDataC = nodeTemp.getDataC();
790     }
791     else
792     nodeDataC = arg.getDataC();
793     #else
794     nodeDataC = arg.getDataC();
795     #endif
796 jgs 82 switch(grad.getFunctionSpace().getTypeCode()) {
797     case(Nodes):
798 jgs 150 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
799 jgs 82 break;
800 gross 1062 case(ReducedNodes):
801     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
802     break;
803 jgs 82 case(Elements):
804 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
805 jgs 82 break;
806 gross 1062 case(ReducedElements):
807     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
808     break;
809 jgs 82 case(FaceElements):
810 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
811 jgs 82 break;
812 gross 1062 case(ReducedFaceElements):
813     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
814     break;
815 jgs 82 case(Points):
816 jgs 150 throw FinleyAdapterException("Error - Gradient at points is not supported.");
817 jgs 82 break;
818     case(ContactElementsZero):
819 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
820 jgs 82 break;
821 gross 1062 case(ReducedContactElementsZero):
822     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
823     break;
824 jgs 82 case(ContactElementsOne):
825 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
826 jgs 82 break;
827 gross 1062 case(ReducedContactElementsOne):
828     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
829     break;
830 jgs 82 case(DegreesOfFreedom):
831 jgs 150 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
832 jgs 82 break;
833     case(ReducedDegreesOfFreedom):
834 jgs 150 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
835 jgs 82 break;
836     default:
837 jgs 150 stringstream temp;
838     temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
839     throw FinleyAdapterException(temp.str());
840 jgs 82 break;
841     }
842     checkFinleyError();
843     }
844 jgs 149
845 jgs 82 //
846     // returns the size of elements:
847     //
848 jgs 480 void MeshAdapter::setToSize(escript::Data& size) const
849 jgs 82 {
850     Finley_Mesh* mesh=m_finleyMesh.get();
851     escriptDataC tmp=size.getDataC();
852     switch(size.getFunctionSpace().getTypeCode()) {
853     case(Nodes):
854 jgs 150 throw FinleyAdapterException("Error - Size of nodes is not supported.");
855 jgs 82 break;
856 gross 1062 case(ReducedNodes):
857     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
858     break;
859 jgs 82 case(Elements):
860     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
861     break;
862 gross 1062 case(ReducedElements):
863     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
864     break;
865 jgs 82 case(FaceElements):
866     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
867     break;
868 gross 1062 case(ReducedFaceElements):
869     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
870     break;
871 jgs 82 case(Points):
872 jgs 150 throw FinleyAdapterException("Error - Size of point elements is not supported.");
873 jgs 82 break;
874     case(ContactElementsZero):
875     case(ContactElementsOne):
876     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
877     break;
878 gross 1062 case(ReducedContactElementsZero):
879     case(ReducedContactElementsOne):
880     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
881     break;
882 jgs 82 case(DegreesOfFreedom):
883 jgs 150 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
884 jgs 82 break;
885     case(ReducedDegreesOfFreedom):
886 jgs 150 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
887 jgs 82 break;
888     default:
889 jgs 150 stringstream temp;
890     temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
891     throw FinleyAdapterException(temp.str());
892 jgs 82 break;
893     }
894     checkFinleyError();
895     }
896 jgs 149
897 jgs 82 // sets the location of nodes:
898 jgs 480 void MeshAdapter::setNewX(const escript::Data& new_x)
899 jgs 82 {
900     Finley_Mesh* mesh=m_finleyMesh.get();
901 bcumming 751 escriptDataC tmp;
902 jgs 82 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
903     if (newDomain!=*this)
904     throw FinleyAdapterException("Error - Illegal domain of new point locations");
905 bcumming 751 tmp = new_x.getDataC();
906     Finley_Mesh_setCoordinates(mesh,&tmp);
907 jgs 82 checkFinleyError();
908     }
909 jgs 149
910 jgs 82 // saves a data array in openDX format:
911 jgs 153 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
912     {
913     int MAX_namelength=256;
914     const int num_data=boost::python::extract<int>(arg.attr("__len__")());
915 woo409 757 /* win32 refactor */
916     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
917     for(int i=0;i<num_data;i++)
918     {
919     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
920     }
921 jgs 153
922 woo409 757 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
923     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
924     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
925    
926 jgs 153 boost::python::list keys=arg.keys();
927     for (int i=0;i<num_data;++i) {
928 jgs 480 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
929 jgs 153 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
930 dhawcroft 793 throw FinleyAdapterException("Error in saveDX: Data must be defined on same Domain");
931 jgs 153 data[i]=d.getDataC();
932     ptr_data[i]=&(data[i]);
933     std::string n=boost::python::extract<std::string>(keys[i]);
934     c_names[i]=&(names[i][0]);
935     if (n.length()>MAX_namelength-1) {
936     strncpy(c_names[i],n.c_str(),MAX_namelength-1);
937     c_names[i][MAX_namelength-1]='\0';
938     } else {
939     strcpy(c_names[i],n.c_str());
940     }
941     }
942     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
943     checkFinleyError();
944 woo409 757
945     /* win32 refactor */
946     TMPMEMFREE(c_names);
947     TMPMEMFREE(data);
948     TMPMEMFREE(ptr_data);
949     for(int i=0;i<num_data;i++)
950     {
951     TMPMEMFREE(names[i]);
952     }
953     TMPMEMFREE(names);
954    
955 jgs 153 return;
956 jgs 82 }
957 jgs 149
958 jgs 110 // saves a data array in openVTK format:
959 jgs 153 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
960     {
961     int MAX_namelength=256;
962     const int num_data=boost::python::extract<int>(arg.attr("__len__")());
963 woo409 757 /* win32 refactor */
964     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
965     for(int i=0;i<num_data;i++)
966     {
967     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
968     }
969 jgs 153
970 woo409 757 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
971     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
972     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
973    
974 jgs 153 boost::python::list keys=arg.keys();
975     for (int i=0;i<num_data;++i) {
976 jgs 480 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
977 jgs 153 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
978     throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
979     data[i]=d.getDataC();
980     ptr_data[i]=&(data[i]);
981     std::string n=boost::python::extract<std::string>(keys[i]);
982     c_names[i]=&(names[i][0]);
983     if (n.length()>MAX_namelength-1) {
984     strncpy(c_names[i],n.c_str(),MAX_namelength-1);
985     c_names[i][MAX_namelength-1]='\0';
986     } else {
987     strcpy(c_names[i],n.c_str());
988     }
989     }
990 dhawcroft 793 #ifndef PASO_MPI
991 jgs 153 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
992 dhawcroft 793 #else
993     Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
994     #endif
995    
996     checkFinleyError();
997 woo409 757 /* win32 refactor */
998     TMPMEMFREE(c_names);
999     TMPMEMFREE(data);
1000     TMPMEMFREE(ptr_data);
1001     for(int i=0;i<num_data;i++)
1002     {
1003     TMPMEMFREE(names[i]);
1004     }
1005     TMPMEMFREE(names);
1006    
1007 jgs 153 return;
1008 jgs 110 }
1009 jgs 153
1010    
1011 jgs 82 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1012     SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1013     const int row_blocksize,
1014     const escript::FunctionSpace& row_functionspace,
1015     const int column_blocksize,
1016     const escript::FunctionSpace& column_functionspace,
1017 jgs 102 const int type) const
1018 jgs 82 {
1019     int reduceRowOrder=0;
1020     int reduceColOrder=0;
1021     // is the domain right?
1022     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1023     if (row_domain!=*this)
1024     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1025     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1026     if (col_domain!=*this)
1027     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1028     // is the function space type right
1029     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1030     reduceRowOrder=0;
1031     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1032     reduceRowOrder=1;
1033     } else {
1034     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1035     }
1036     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1037     reduceColOrder=0;
1038     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1039     reduceColOrder=1;
1040     } else {
1041     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1042     }
1043     // generate matrix:
1044 jgs 102
1045 jgs 150 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1046 jgs 82 checkFinleyError();
1047 ksteube 971 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1048 jgs 150 checkPasoError();
1049     Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
1050 jgs 82 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1051     }
1052 jgs 149
1053 jgs 82 //
1054     // vtkObject MeshAdapter::createVtkObject() const
1055     // TODO:
1056     //
1057 jgs 149
1058 jgs 82 //
1059     // returns true if data at the atom_type is considered as being cell centered:
1060     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1061     {
1062     switch(functionSpaceCode) {
1063     case(Nodes):
1064     case(DegreesOfFreedom):
1065     case(ReducedDegreesOfFreedom):
1066     return false;
1067     break;
1068     case(Elements):
1069     case(FaceElements):
1070     case(Points):
1071     case(ContactElementsZero):
1072     case(ContactElementsOne):
1073 gross 1062 case(ReducedElements):
1074     case(ReducedFaceElements):
1075     case(ReducedContactElementsZero):
1076     case(ReducedContactElementsOne):
1077 jgs 82 return true;
1078     break;
1079     default:
1080 jgs 150 stringstream temp;
1081     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1082     throw FinleyAdapterException(temp.str());
1083 jgs 82 break;
1084     }
1085     checkFinleyError();
1086     return false;
1087     }
1088 jgs 149
1089 jgs 82 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1090     {
1091     switch(functionSpaceType_source) {
1092     case(Nodes):
1093     switch(functionSpaceType_target) {
1094     case(Nodes):
1095 gross 1062 case(ReducedNodes):
1096 jgs 82 case(ReducedDegreesOfFreedom):
1097     case(DegreesOfFreedom):
1098     case(Elements):
1099 gross 1062 case(ReducedElements):
1100 jgs 82 case(FaceElements):
1101 gross 1062 case(ReducedFaceElements):
1102 jgs 82 case(Points):
1103     case(ContactElementsZero):
1104 gross 1062 case(ReducedContactElementsZero):
1105 jgs 82 case(ContactElementsOne):
1106 gross 1062 case(ReducedContactElementsOne):
1107 jgs 82 return true;
1108     default:
1109 jgs 150 stringstream temp;
1110     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1111     throw FinleyAdapterException(temp.str());
1112 jgs 82 }
1113     break;
1114 gross 1062 case(ReducedNodes):
1115     switch(functionSpaceType_target) {
1116     case(ReducedNodes):
1117     case(ReducedDegreesOfFreedom):
1118     case(Elements):
1119     case(ReducedElements):
1120     case(FaceElements):
1121     case(ReducedFaceElements):
1122     case(Points):
1123     case(ContactElementsZero):
1124     case(ReducedContactElementsZero):
1125     case(ContactElementsOne):
1126     case(ReducedContactElementsOne):
1127     return true;
1128     case(Nodes):
1129     case(DegreesOfFreedom):
1130     return false;
1131     default:
1132     stringstream temp;
1133     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1134     throw FinleyAdapterException(temp.str());
1135     }
1136     break;
1137 jgs 82 case(Elements):
1138     if (functionSpaceType_target==Elements) {
1139     return true;
1140 gross 1116 } else if (functionSpaceType_target==ReducedElements) {
1141     return true;
1142 jgs 82 } else {
1143     return false;
1144     }
1145 gross 1062 case(ReducedElements):
1146     if (functionSpaceType_target==ReducedElements) {
1147     return true;
1148     } else {
1149     return false;
1150     }
1151 jgs 82 case(FaceElements):
1152     if (functionSpaceType_target==FaceElements) {
1153     return true;
1154 gross 1116 } else if (functionSpaceType_target==ReducedFaceElements) {
1155     return true;
1156 jgs 82 } else {
1157     return false;
1158     }
1159 gross 1062 case(ReducedFaceElements):
1160     if (functionSpaceType_target==ReducedFaceElements) {
1161     return true;
1162     } else {
1163     return false;
1164     }
1165 jgs 82 case(Points):
1166     if (functionSpaceType_target==Points) {
1167     return true;
1168     } else {
1169     return false;
1170     }
1171     case(ContactElementsZero):
1172     case(ContactElementsOne):
1173     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1174     return true;
1175 gross 1116 } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1176     return true;
1177 jgs 82 } else {
1178     return false;
1179     }
1180 gross 1062 case(ReducedContactElementsZero):
1181     case(ReducedContactElementsOne):
1182     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1183     return true;
1184     } else {
1185     return false;
1186     }
1187 jgs 82 case(DegreesOfFreedom):
1188     switch(functionSpaceType_target) {
1189     case(ReducedDegreesOfFreedom):
1190     case(DegreesOfFreedom):
1191     case(Nodes):
1192 gross 1062 case(ReducedNodes):
1193 jgs 82 case(Elements):
1194 gross 1062 case(ReducedElements):
1195     case(Points):
1196 jgs 82 case(FaceElements):
1197 gross 1062 case(ReducedFaceElements):
1198 jgs 82 case(ContactElementsZero):
1199 gross 1062 case(ReducedContactElementsZero):
1200 jgs 82 case(ContactElementsOne):
1201 gross 1062 case(ReducedContactElementsOne):
1202 jgs 82 return true;
1203     default:
1204 jgs 150 stringstream temp;
1205     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1206     throw FinleyAdapterException(temp.str());
1207 jgs 82 }
1208     break;
1209     case(ReducedDegreesOfFreedom):
1210     switch(functionSpaceType_target) {
1211     case(ReducedDegreesOfFreedom):
1212 gross 1062 case(ReducedNodes):
1213 jgs 82 case(Elements):
1214 gross 1062 case(ReducedElements):
1215 jgs 82 case(FaceElements):
1216 gross 1062 case(ReducedFaceElements):
1217 jgs 82 case(Points):
1218     case(ContactElementsZero):
1219 gross 1062 case(ReducedContactElementsZero):
1220 jgs 82 case(ContactElementsOne):
1221 gross 1062 case(ReducedContactElementsOne):
1222 jgs 82 return true;
1223 jgs 147 case(Nodes):
1224 jgs 82 case(DegreesOfFreedom):
1225     return false;
1226     default:
1227 jgs 150 stringstream temp;
1228     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1229     throw FinleyAdapterException(temp.str());
1230 jgs 82 }
1231     break;
1232     default:
1233 jgs 150 stringstream temp;
1234     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1235     throw FinleyAdapterException(temp.str());
1236 jgs 82 break;
1237     }
1238     checkFinleyError();
1239     return false;
1240     }
1241 jgs 149
1242 jgs 82 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1243     {
1244     return false;
1245     }
1246 jgs 104
1247 jgs 121 bool MeshAdapter::operator==(const AbstractDomain& other) const
1248 jgs 82 {
1249 jgs 121 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1250     if (temp!=0) {
1251     return (m_finleyMesh==temp->m_finleyMesh);
1252     } else {
1253     return false;
1254     }
1255 jgs 82 }
1256    
1257 jgs 121 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1258 jgs 104 {
1259 jgs 121 return !(operator==(other));
1260 jgs 104 }
1261    
1262 jgs 150 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1263 jgs 102 {
1264 jgs 150 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1265     checkPasoError();
1266 jgs 102 return out;
1267     }
1268 jgs 149
1269 jgs 480 escript::Data MeshAdapter::getX() const
1270 jgs 102 {
1271     return continuousFunction(asAbstractContinuousDomain()).getX();
1272     }
1273 jgs 149
1274 jgs 480 escript::Data MeshAdapter::getNormal() const
1275 jgs 102 {
1276     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1277     }
1278 jgs 149
1279 jgs 480 escript::Data MeshAdapter::getSize() const
1280 jgs 102 {
1281     return function(asAbstractContinuousDomain()).getSize();
1282     }
1283    
1284 gross 1062 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1285     {
1286     int *out=0,i;
1287     Finley_Mesh* mesh=m_finleyMesh.get();
1288     switch (functionSpaceType) {
1289     case(Nodes):
1290     if (mesh->Nodes!=NULL) {
1291     out=mesh->Nodes->Id;
1292     break;
1293     }
1294     case(ReducedNodes):
1295     throw FinleyAdapterException("Error - ReducedNodes not supported yet.");
1296     break;
1297     case(Elements):
1298 gross 1080 out=mesh->Elements->Id;
1299 gross 1062 break;
1300     case(ReducedElements):
1301     out=mesh->Elements->Id;
1302     break;
1303     case(FaceElements):
1304     out=mesh->FaceElements->Id;
1305     break;
1306     case(ReducedFaceElements):
1307     out=mesh->FaceElements->Id;
1308     break;
1309     case(Points):
1310     out=mesh->Points->Id;
1311     break;
1312     case(ContactElementsZero):
1313     out=mesh->ContactElements->Id;
1314     break;
1315     case(ReducedContactElementsZero):
1316     out=mesh->ContactElements->Id;
1317     break;
1318     case(ContactElementsOne):
1319     out=mesh->ContactElements->Id;
1320     break;
1321     case(ReducedContactElementsOne):
1322     out=mesh->ContactElements->Id;
1323     break;
1324     case(DegreesOfFreedom):
1325     out=mesh->Nodes->degreeOfFreedomId;
1326     break;
1327     case(ReducedDegreesOfFreedom):
1328     out=mesh->Nodes->reducedDegreeOfFreedomId;
1329     break;
1330     default:
1331     stringstream temp;
1332     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1333     throw FinleyAdapterException(temp.str());
1334     break;
1335     }
1336     return out;
1337     }
1338 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1339     {
1340 gross 547 int out=0;
1341     Finley_Mesh* mesh=m_finleyMesh.get();
1342     switch (functionSpaceType) {
1343     case(Nodes):
1344     out=mesh->Nodes->Tag[sampleNo];
1345     break;
1346 gross 1062 case(ReducedNodes):
1347     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1348     break;
1349 gross 547 case(Elements):
1350     out=mesh->Elements->Tag[sampleNo];
1351     break;
1352 gross 1062 case(ReducedElements):
1353     out=mesh->Elements->Tag[sampleNo];
1354     break;
1355 gross 547 case(FaceElements):
1356     out=mesh->FaceElements->Tag[sampleNo];
1357     break;
1358 gross 1062 case(ReducedFaceElements):
1359     out=mesh->FaceElements->Tag[sampleNo];
1360     break;
1361 gross 547 case(Points):
1362     out=mesh->Points->Tag[sampleNo];
1363     break;
1364     case(ContactElementsZero):
1365     out=mesh->ContactElements->Tag[sampleNo];
1366     break;
1367 gross 1062 case(ReducedContactElementsZero):
1368     out=mesh->ContactElements->Tag[sampleNo];
1369     break;
1370 gross 547 case(ContactElementsOne):
1371     out=mesh->ContactElements->Tag[sampleNo];
1372     break;
1373 gross 1062 case(ReducedContactElementsOne):
1374     out=mesh->ContactElements->Tag[sampleNo];
1375 gross 547 break;
1376     case(DegreesOfFreedom):
1377 gross 1062 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1378 gross 547 break;
1379     case(ReducedDegreesOfFreedom):
1380 gross 1062 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1381 gross 547 break;
1382     default:
1383     stringstream temp;
1384     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1385     throw FinleyAdapterException(temp.str());
1386     break;
1387     }
1388     return out;
1389 jgs 110 }
1390    
1391 gross 1062
1392 gross 767 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1393     {
1394     Finley_Mesh* mesh=m_finleyMesh.get();
1395     escriptDataC tmp=mask.getDataC();
1396     switch(functionSpaceType) {
1397     case(Nodes):
1398     Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1399     break;
1400 gross 1062 case(ReducedNodes):
1401     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1402     break;
1403 gross 767 case(DegreesOfFreedom):
1404     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1405     break;
1406     case(ReducedDegreesOfFreedom):
1407     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1408     break;
1409     case(Elements):
1410     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1411     break;
1412 gross 1062 case(ReducedElements):
1413     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1414     break;
1415 gross 767 case(FaceElements):
1416     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1417     break;
1418 gross 1062 case(ReducedFaceElements):
1419     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1420     break;
1421 gross 767 case(Points):
1422     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1423     break;
1424     case(ContactElementsZero):
1425     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1426     break;
1427 gross 1062 case(ReducedContactElementsZero):
1428     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1429     break;
1430 gross 767 case(ContactElementsOne):
1431     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1432     break;
1433 gross 1062 case(ReducedContactElementsOne):
1434     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1435     break;
1436 gross 767 default:
1437     stringstream temp;
1438     temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1439     throw FinleyAdapterException(temp.str());
1440     }
1441     checkFinleyError();
1442     return;
1443     }
1444    
1445 gross 1044 void MeshAdapter::setTagMap(const std::string& name, int tag)
1446     {
1447     Finley_Mesh* mesh=m_finleyMesh.get();
1448     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1449     checkPasoError();
1450     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1451     }
1452 gross 767
1453 gross 1044 int MeshAdapter::getTag(const std::string& name) const
1454     {
1455     Finley_Mesh* mesh=m_finleyMesh.get();
1456     int tag=0;
1457     tag=Finley_Mesh_getTag(mesh, name.c_str());
1458     checkPasoError();
1459     // throwStandardException("MeshAdapter::getTag is not implemented.");
1460     return tag;
1461     }
1462    
1463     bool MeshAdapter::isValidTagName(const std::string& name) const
1464     {
1465     Finley_Mesh* mesh=m_finleyMesh.get();
1466     return Finley_Mesh_isValidTagName(mesh,name.c_str());
1467     }
1468    
1469     std::string MeshAdapter::showTagNames() const
1470     {
1471     stringstream temp;
1472     Finley_Mesh* mesh=m_finleyMesh.get();
1473     Finley_TagMap* tag_map=mesh->TagMap;
1474     while (tag_map) {
1475     temp << tag_map->name;
1476     tag_map=tag_map->next;
1477     if (tag_map) temp << ", ";
1478     }
1479 gross 1137 return string(temp.str());
1480 gross 1044 }
1481    
1482 jgs 82 } // end of namespace

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26