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

Annotation of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1062 - (hide annotations)
Mon Mar 26 06:17:53 2007 UTC (12 years ago) by gross
File size: 55217 byte(s)
reduced integration schemes are implemented now for grad, integrate, etc. Tests still to be added.
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     return "FinleyMesh";
91     }
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     return loc->second;
101     }
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 jgs 82 } else {
458 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
459 jgs 82 }
460     break;
461 gross 1062 case(ReducedElements):
462     if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
463     Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
464     } else {
465     throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
466     }
467     break;
468 jgs 82 case(FaceElements):
469     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
470 gross 1062 Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
471 jgs 82 } else {
472 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
473 jgs 82 }
474 gross 1062 break;
475     case(ReducedFaceElements):
476     if (target.getFunctionSpace().getTypeCode()==FaceElements) {
477     Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
478     } else {
479     throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
480     }
481     break;
482 jgs 82 case(Points):
483     if (target.getFunctionSpace().getTypeCode()==Points) {
484 gross 1062 Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
485 jgs 82 } else {
486 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on points possible.");
487 jgs 82 }
488     break;
489     case(ContactElementsZero):
490     case(ContactElementsOne):
491     if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
492 gross 1062 Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
493 jgs 82 } else {
494 jgs 150 throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
495 jgs 82 break;
496     }
497     break;
498 gross 1062 case(ReducedContactElementsZero):
499     case(ReducedContactElementsOne):
500     if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
501     Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
502     } else {
503     throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
504     break;
505     }
506     break;
507 bcumming 751 case(DegreesOfFreedom):
508 jgs 82 switch(target.getFunctionSpace().getTypeCode()) {
509     case(ReducedDegreesOfFreedom):
510     case(DegreesOfFreedom):
511     case(Nodes):
512 gross 1062 case(ReducedNodes):
513     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
514 jgs 82 break;
515 bcumming 751 #ifndef PASO_MPI
516 jgs 82 case(Elements):
517 gross 1062 case(ReducedElements):
518     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
519 jgs 82 break;
520     case(FaceElements):
521 gross 1062 case(ReducedFaceElements):
522     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
523 jgs 82 break;
524     case(Points):
525 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
526 jgs 82 break;
527     case(ContactElementsZero):
528     case(ContactElementsOne):
529 gross 1062 case(ReducedContactElementsZero):
530     case(ReducedContactElementsOne):
531     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
532 jgs 82 break;
533 bcumming 751 #else
534     /* need to copy Degrees of freedom data to nodal data so that the external values are available */
535     case(Elements):
536 gross 1062 case(ReducedElements):
537 bcumming 751 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
538 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);
539 bcumming 751 break;
540     case(FaceElements):
541 gross 1062 case(ReducedFaceElements):
542 bcumming 751 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
543 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);
544 bcumming 751 break;
545     case(Points):
546     escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
547 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);
548 bcumming 751 break;
549     case(ContactElementsZero):
550     case(ContactElementsOne):
551 gross 1062 case(ReducedContactElementsZero):
552     case(ReducedContactElementsOne):
553 bcumming 751 escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
554 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);
555 bcumming 751 break;
556     #endif
557 jgs 82 default:
558 jgs 150 stringstream temp;
559     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
560     throw FinleyAdapterException(temp.str());
561 jgs 82 break;
562     }
563     break;
564     case(ReducedDegreesOfFreedom):
565     switch(target.getFunctionSpace().getTypeCode()) {
566 gross 1062 case(Nodes):
567     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
568     break;
569     case(ReducedNodes):
570     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
571     break;
572     case(DegreesOfFreedom):
573     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
574     break;
575 jgs 82 case(ReducedDegreesOfFreedom):
576 gross 1062 Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
577 jgs 82 break;
578     case(Elements):
579 gross 1062 case(ReducedElements):
580     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
581 jgs 82 break;
582     case(FaceElements):
583 gross 1062 case(ReducedFaceElements):
584     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
585 jgs 82 break;
586     case(Points):
587 gross 1062 Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
588 jgs 82 break;
589     case(ContactElementsZero):
590     case(ContactElementsOne):
591 gross 1062 case(ReducedContactElementsZero):
592     case(ReducedContactElementsOne):
593     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
594 jgs 82 break;
595     default:
596 jgs 150 stringstream temp;
597     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
598     throw FinleyAdapterException(temp.str());
599 jgs 82 break;
600     }
601     break;
602     default:
603 jgs 150 stringstream temp;
604     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
605     throw FinleyAdapterException(temp.str());
606 jgs 82 break;
607     }
608     checkFinleyError();
609     }
610    
611     //
612     // copies the locations of sample points into x:
613     //
614 jgs 480 void MeshAdapter::setToX(escript::Data& arg) const
615 jgs 82 {
616     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
617     if (argDomain!=*this)
618     throw FinleyAdapterException("Error - Illegal domain of data point locations");
619     Finley_Mesh* mesh=m_finleyMesh.get();
620     // in case of values node coordinates we can do the job directly:
621     if (arg.getFunctionSpace().getTypeCode()==Nodes) {
622 gross 1062 escriptDataC _arg=arg.getDataC();
623     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
624 jgs 82 } else {
625 jgs 480 escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);
626 gross 1062 escriptDataC _tmp_data=tmp_data.getDataC();
627     Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
628 jgs 82 // this is then interpolated onto arg:
629     interpolateOnDomain(arg,tmp_data);
630     }
631     checkFinleyError();
632     }
633 jgs 149
634 jgs 82 //
635     // return the normal vectors at the location of data points as a Data object:
636     //
637 jgs 480 void MeshAdapter::setToNormal(escript::Data& normal) const
638 jgs 82 {
639     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());
640     if (normalDomain!=*this)
641     throw FinleyAdapterException("Error - Illegal domain of normal locations");
642     Finley_Mesh* mesh=m_finleyMesh.get();
643 gross 1062 escriptDataC _normal=normal.getDataC();
644 jgs 82 switch(normal.getFunctionSpace().getTypeCode()) {
645     case(Nodes):
646 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
647 jgs 82 break;
648 gross 1062 case(ReducedNodes):
649     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced nodes");
650     break;
651 jgs 82 case(Elements):
652 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements");
653 jgs 82 break;
654 gross 1062 case(ReducedElements):
655     throw FinleyAdapterException("Error - Finley does not support surface normal vectors for elements with reduced integration order");
656     break;
657 jgs 82 case (FaceElements):
658 gross 1062 Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
659 jgs 82 break;
660 gross 1062 case (ReducedFaceElements):
661     Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
662     break;
663 jgs 82 case(Points):
664 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
665 jgs 82 break;
666     case (ContactElementsOne):
667     case (ContactElementsZero):
668 gross 1062 Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
669 jgs 82 break;
670 gross 1062 case (ReducedContactElementsOne):
671     case (ReducedContactElementsZero):
672     Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);
673     break;
674 jgs 82 case(DegreesOfFreedom):
675 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
676 jgs 82 break;
677     case(ReducedDegreesOfFreedom):
678 jgs 150 throw FinleyAdapterException("Error - Finley does not support surface normal vectors for reduced degrees of freedom.");
679 jgs 82 break;
680     default:
681 jgs 150 stringstream temp;
682     temp << "Error - Normal Vectors: Finley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
683     throw FinleyAdapterException(temp.str());
684 jgs 82 break;
685     }
686     checkFinleyError();
687     }
688 jgs 149
689 jgs 82 //
690     // interpolates data to other domain:
691     //
692 jgs 480 void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
693 jgs 82 {
694     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
695     if (targetDomain!=*this)
696     throw FinleyAdapterException("Error - Illegal domain of interpolation target");
697    
698 jgs 150 throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
699 jgs 82 }
700 jgs 149
701 jgs 82 //
702     // calculates the integral of a function defined of arg:
703     //
704 jgs 480 void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
705 jgs 82 {
706     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
707     if (argDomain!=*this)
708     throw FinleyAdapterException("Error - Illegal domain of integration kernel");
709    
710     Finley_Mesh* mesh=m_finleyMesh.get();
711 gross 1062 escriptDataC _arg=arg.getDataC();
712 jgs 82 switch(arg.getFunctionSpace().getTypeCode()) {
713     case(Nodes):
714 gross 1062 /* TODO */
715 jgs 150 throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");
716 jgs 82 break;
717 gross 1062 case(ReducedNodes):
718     /* TODO */
719     throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");
720     break;
721 jgs 82 case(Elements):
722 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
723 jgs 82 break;
724 gross 1062 case(ReducedElements):
725     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
726 jgs 82 case(FaceElements):
727 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
728 jgs 82 break;
729 gross 1062 case(ReducedFaceElements):
730     Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
731     break;
732 jgs 82 case(Points):
733 jgs 150 throw FinleyAdapterException("Error - Integral of data on points is not supported.");
734 jgs 82 break;
735     case(ContactElementsZero):
736 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
737 jgs 82 break;
738 gross 1062 case(ReducedContactElementsZero):
739     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
740     break;
741 jgs 82 case(ContactElementsOne):
742 gross 1062 Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
743 jgs 82 break;
744 gross 1062 case(ReducedContactElementsOne):
745     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
746     break;
747 jgs 82 case(DegreesOfFreedom):
748 jgs 150 throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");
749 jgs 82 break;
750     case(ReducedDegreesOfFreedom):
751 jgs 150 throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");
752 jgs 82 break;
753     default:
754 jgs 150 stringstream temp;
755     temp << "Error - Integrals: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
756     throw FinleyAdapterException(temp.str());
757 jgs 82 break;
758     }
759     checkFinleyError();
760     }
761 jgs 149
762 jgs 82 //
763     // calculates the gradient of arg:
764     //
765 jgs 480 void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
766 jgs 82 {
767     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());
768     if (argDomain!=*this)
769     throw FinleyAdapterException("Error - Illegal domain of gradient argument");
770     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());
771     if (gradDomain!=*this)
772     throw FinleyAdapterException("Error - Illegal domain of gradient");
773    
774     Finley_Mesh* mesh=m_finleyMesh.get();
775 gross 1062 escriptDataC _grad=grad.getDataC();
776 bcumming 751 escriptDataC nodeDataC;
777     #ifdef PASO_MPI
778     escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
779     if( arg.getFunctionSpace().getTypeCode() != Nodes )
780     {
781     Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
782     nodeDataC = nodeTemp.getDataC();
783     }
784     else
785     nodeDataC = arg.getDataC();
786     #else
787     nodeDataC = arg.getDataC();
788     #endif
789 jgs 82 switch(grad.getFunctionSpace().getTypeCode()) {
790     case(Nodes):
791 jgs 150 throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
792 jgs 82 break;
793 gross 1062 case(ReducedNodes):
794     throw FinleyAdapterException("Error - Gradient at reduced nodes is not supported.");
795     break;
796 jgs 82 case(Elements):
797 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
798 jgs 82 break;
799 gross 1062 case(ReducedElements):
800     Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
801     break;
802 jgs 82 case(FaceElements):
803 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
804 jgs 82 break;
805 gross 1062 case(ReducedFaceElements):
806     Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
807     break;
808 jgs 82 case(Points):
809 jgs 150 throw FinleyAdapterException("Error - Gradient at points is not supported.");
810 jgs 82 break;
811     case(ContactElementsZero):
812 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
813 jgs 82 break;
814 gross 1062 case(ReducedContactElementsZero):
815     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
816     break;
817 jgs 82 case(ContactElementsOne):
818 gross 1062 Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
819 jgs 82 break;
820 gross 1062 case(ReducedContactElementsOne):
821     Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);
822     break;
823 jgs 82 case(DegreesOfFreedom):
824 jgs 150 throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
825 jgs 82 break;
826     case(ReducedDegreesOfFreedom):
827 jgs 150 throw FinleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
828 jgs 82 break;
829     default:
830 jgs 150 stringstream temp;
831     temp << "Error - Gradient: Finley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
832     throw FinleyAdapterException(temp.str());
833 jgs 82 break;
834     }
835     checkFinleyError();
836     }
837 jgs 149
838 jgs 82 //
839     // returns the size of elements:
840     //
841 jgs 480 void MeshAdapter::setToSize(escript::Data& size) const
842 jgs 82 {
843     Finley_Mesh* mesh=m_finleyMesh.get();
844     escriptDataC tmp=size.getDataC();
845     switch(size.getFunctionSpace().getTypeCode()) {
846     case(Nodes):
847 jgs 150 throw FinleyAdapterException("Error - Size of nodes is not supported.");
848 jgs 82 break;
849 gross 1062 case(ReducedNodes):
850     throw FinleyAdapterException("Error - Size of reduced nodes is not supported.");
851     break;
852 jgs 82 case(Elements):
853     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
854     break;
855 gross 1062 case(ReducedElements):
856     Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
857     break;
858 jgs 82 case(FaceElements):
859     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
860     break;
861 gross 1062 case(ReducedFaceElements):
862     Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
863     break;
864 jgs 82 case(Points):
865 jgs 150 throw FinleyAdapterException("Error - Size of point elements is not supported.");
866 jgs 82 break;
867     case(ContactElementsZero):
868     case(ContactElementsOne):
869     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
870     break;
871 gross 1062 case(ReducedContactElementsZero):
872     case(ReducedContactElementsOne):
873     Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);
874     break;
875 jgs 82 case(DegreesOfFreedom):
876 jgs 150 throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
877 jgs 82 break;
878     case(ReducedDegreesOfFreedom):
879 jgs 150 throw FinleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
880 jgs 82 break;
881     default:
882 jgs 150 stringstream temp;
883     temp << "Error - Element size: Finley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
884     throw FinleyAdapterException(temp.str());
885 jgs 82 break;
886     }
887     checkFinleyError();
888     }
889 jgs 149
890 jgs 82 // sets the location of nodes:
891 jgs 480 void MeshAdapter::setNewX(const escript::Data& new_x)
892 jgs 82 {
893     Finley_Mesh* mesh=m_finleyMesh.get();
894 bcumming 751 escriptDataC tmp;
895 jgs 82 const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
896     if (newDomain!=*this)
897     throw FinleyAdapterException("Error - Illegal domain of new point locations");
898 bcumming 751 tmp = new_x.getDataC();
899     Finley_Mesh_setCoordinates(mesh,&tmp);
900 jgs 82 checkFinleyError();
901     }
902 jgs 149
903 jgs 82 // saves a data array in openDX format:
904 jgs 153 void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
905     {
906     int MAX_namelength=256;
907     const int num_data=boost::python::extract<int>(arg.attr("__len__")());
908 woo409 757 /* win32 refactor */
909     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
910     for(int i=0;i<num_data;i++)
911     {
912     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
913     }
914 jgs 153
915 woo409 757 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
916     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
917     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
918    
919 jgs 153 boost::python::list keys=arg.keys();
920     for (int i=0;i<num_data;++i) {
921 jgs 480 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
922 jgs 153 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
923 dhawcroft 793 throw FinleyAdapterException("Error in saveDX: Data must be defined on same Domain");
924 jgs 153 data[i]=d.getDataC();
925     ptr_data[i]=&(data[i]);
926     std::string n=boost::python::extract<std::string>(keys[i]);
927     c_names[i]=&(names[i][0]);
928     if (n.length()>MAX_namelength-1) {
929     strncpy(c_names[i],n.c_str(),MAX_namelength-1);
930     c_names[i][MAX_namelength-1]='\0';
931     } else {
932     strcpy(c_names[i],n.c_str());
933     }
934     }
935     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
936     checkFinleyError();
937 woo409 757
938     /* win32 refactor */
939     TMPMEMFREE(c_names);
940     TMPMEMFREE(data);
941     TMPMEMFREE(ptr_data);
942     for(int i=0;i<num_data;i++)
943     {
944     TMPMEMFREE(names[i]);
945     }
946     TMPMEMFREE(names);
947    
948 jgs 153 return;
949 jgs 82 }
950 jgs 149
951 jgs 110 // saves a data array in openVTK format:
952 jgs 153 void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
953     {
954     int MAX_namelength=256;
955     const int num_data=boost::python::extract<int>(arg.attr("__len__")());
956 woo409 757 /* win32 refactor */
957     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
958     for(int i=0;i<num_data;i++)
959     {
960     names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;
961     }
962 jgs 153
963 woo409 757 char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
964     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
965     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
966    
967 jgs 153 boost::python::list keys=arg.keys();
968     for (int i=0;i<num_data;++i) {
969 jgs 480 escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
970 jgs 153 if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)
971     throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
972     data[i]=d.getDataC();
973     ptr_data[i]=&(data[i]);
974     std::string n=boost::python::extract<std::string>(keys[i]);
975     c_names[i]=&(names[i][0]);
976     if (n.length()>MAX_namelength-1) {
977     strncpy(c_names[i],n.c_str(),MAX_namelength-1);
978     c_names[i][MAX_namelength-1]='\0';
979     } else {
980     strcpy(c_names[i],n.c_str());
981     }
982     }
983 dhawcroft 793 #ifndef PASO_MPI
984 jgs 153 Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
985 dhawcroft 793 #else
986     Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
987     #endif
988    
989     checkFinleyError();
990 woo409 757 /* win32 refactor */
991     TMPMEMFREE(c_names);
992     TMPMEMFREE(data);
993     TMPMEMFREE(ptr_data);
994     for(int i=0;i<num_data;i++)
995     {
996     TMPMEMFREE(names[i]);
997     }
998     TMPMEMFREE(names);
999    
1000 jgs 153 return;
1001 jgs 110 }
1002 jgs 153
1003    
1004 jgs 82 // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:
1005     SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1006     const int row_blocksize,
1007     const escript::FunctionSpace& row_functionspace,
1008     const int column_blocksize,
1009     const escript::FunctionSpace& column_functionspace,
1010 jgs 102 const int type) const
1011 jgs 82 {
1012     int reduceRowOrder=0;
1013     int reduceColOrder=0;
1014     // is the domain right?
1015     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());
1016     if (row_domain!=*this)
1017     throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1018     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());
1019     if (col_domain!=*this)
1020     throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1021     // is the function space type right
1022     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1023     reduceRowOrder=0;
1024     } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1025     reduceRowOrder=1;
1026     } else {
1027     throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1028     }
1029     if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
1030     reduceColOrder=0;
1031     } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1032     reduceColOrder=1;
1033     } else {
1034     throw FinleyAdapterException("Error - illegal function space type for system matrix columns.");
1035     }
1036     // generate matrix:
1037 jgs 102
1038 jgs 150 Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder);
1039 jgs 82 checkFinleyError();
1040 ksteube 971 Paso_SystemMatrix* fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1041 jgs 150 checkPasoError();
1042     Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);
1043 jgs 82 return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1044     }
1045 jgs 149
1046 jgs 82 //
1047     // vtkObject MeshAdapter::createVtkObject() const
1048     // TODO:
1049     //
1050 jgs 149
1051 jgs 82 //
1052     // returns true if data at the atom_type is considered as being cell centered:
1053     bool MeshAdapter::isCellOriented(int functionSpaceCode) const
1054     {
1055     switch(functionSpaceCode) {
1056     case(Nodes):
1057     case(DegreesOfFreedom):
1058     case(ReducedDegreesOfFreedom):
1059     return false;
1060     break;
1061     case(Elements):
1062     case(FaceElements):
1063     case(Points):
1064     case(ContactElementsZero):
1065     case(ContactElementsOne):
1066 gross 1062 case(ReducedElements):
1067     case(ReducedFaceElements):
1068     case(ReducedContactElementsZero):
1069     case(ReducedContactElementsOne):
1070 jgs 82 return true;
1071     break;
1072     default:
1073 jgs 150 stringstream temp;
1074     temp << "Error - Cell: Finley does not know anything about function space type " << functionSpaceCode;
1075     throw FinleyAdapterException(temp.str());
1076 jgs 82 break;
1077     }
1078     checkFinleyError();
1079     return false;
1080     }
1081 jgs 149
1082 jgs 82 bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1083     {
1084     switch(functionSpaceType_source) {
1085     case(Nodes):
1086     switch(functionSpaceType_target) {
1087     case(Nodes):
1088 gross 1062 case(ReducedNodes):
1089 jgs 82 case(ReducedDegreesOfFreedom):
1090     case(DegreesOfFreedom):
1091     case(Elements):
1092 gross 1062 case(ReducedElements):
1093 jgs 82 case(FaceElements):
1094 gross 1062 case(ReducedFaceElements):
1095 jgs 82 case(Points):
1096     case(ContactElementsZero):
1097 gross 1062 case(ReducedContactElementsZero):
1098 jgs 82 case(ContactElementsOne):
1099 gross 1062 case(ReducedContactElementsOne):
1100 jgs 82 return true;
1101     default:
1102 jgs 150 stringstream temp;
1103     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1104     throw FinleyAdapterException(temp.str());
1105 jgs 82 }
1106     break;
1107 gross 1062 case(ReducedNodes):
1108     switch(functionSpaceType_target) {
1109     case(ReducedNodes):
1110     case(ReducedDegreesOfFreedom):
1111     case(Elements):
1112     case(ReducedElements):
1113     case(FaceElements):
1114     case(ReducedFaceElements):
1115     case(Points):
1116     case(ContactElementsZero):
1117     case(ReducedContactElementsZero):
1118     case(ContactElementsOne):
1119     case(ReducedContactElementsOne):
1120     return true;
1121     case(Nodes):
1122     case(DegreesOfFreedom):
1123     return false;
1124     default:
1125     stringstream temp;
1126     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1127     throw FinleyAdapterException(temp.str());
1128     }
1129     break;
1130 jgs 82 case(Elements):
1131     if (functionSpaceType_target==Elements) {
1132     return true;
1133     } else {
1134     return false;
1135     }
1136 gross 1062 case(ReducedElements):
1137     if (functionSpaceType_target==ReducedElements) {
1138     return true;
1139     } else {
1140     return false;
1141     }
1142 jgs 82 case(FaceElements):
1143     if (functionSpaceType_target==FaceElements) {
1144     return true;
1145     } else {
1146     return false;
1147     }
1148 gross 1062 case(ReducedFaceElements):
1149     if (functionSpaceType_target==ReducedFaceElements) {
1150     return true;
1151     } else {
1152     return false;
1153     }
1154 jgs 82 case(Points):
1155     if (functionSpaceType_target==Points) {
1156     return true;
1157     } else {
1158     return false;
1159     }
1160     case(ContactElementsZero):
1161     case(ContactElementsOne):
1162     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1163     return true;
1164     } else {
1165     return false;
1166     }
1167 gross 1062 case(ReducedContactElementsZero):
1168     case(ReducedContactElementsOne):
1169     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1170     return true;
1171     } else {
1172     return false;
1173     }
1174 jgs 82 case(DegreesOfFreedom):
1175     switch(functionSpaceType_target) {
1176     case(ReducedDegreesOfFreedom):
1177     case(DegreesOfFreedom):
1178     case(Nodes):
1179 gross 1062 case(ReducedNodes):
1180 jgs 82 case(Elements):
1181 gross 1062 case(ReducedElements):
1182     case(Points):
1183 jgs 82 case(FaceElements):
1184 gross 1062 case(ReducedFaceElements):
1185 jgs 82 case(ContactElementsZero):
1186 gross 1062 case(ReducedContactElementsZero):
1187 jgs 82 case(ContactElementsOne):
1188 gross 1062 case(ReducedContactElementsOne):
1189 jgs 82 return true;
1190     default:
1191 jgs 150 stringstream temp;
1192     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1193     throw FinleyAdapterException(temp.str());
1194 jgs 82 }
1195     break;
1196     case(ReducedDegreesOfFreedom):
1197     switch(functionSpaceType_target) {
1198     case(ReducedDegreesOfFreedom):
1199 gross 1062 case(ReducedNodes):
1200 jgs 82 case(Elements):
1201 gross 1062 case(ReducedElements):
1202 jgs 82 case(FaceElements):
1203 gross 1062 case(ReducedFaceElements):
1204 jgs 82 case(Points):
1205     case(ContactElementsZero):
1206 gross 1062 case(ReducedContactElementsZero):
1207 jgs 82 case(ContactElementsOne):
1208 gross 1062 case(ReducedContactElementsOne):
1209 jgs 82 return true;
1210 jgs 147 case(Nodes):
1211 jgs 82 case(DegreesOfFreedom):
1212     return false;
1213     default:
1214 jgs 150 stringstream temp;
1215     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1216     throw FinleyAdapterException(temp.str());
1217 jgs 82 }
1218     break;
1219     default:
1220 jgs 150 stringstream temp;
1221     temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
1222     throw FinleyAdapterException(temp.str());
1223 jgs 82 break;
1224     }
1225     checkFinleyError();
1226     return false;
1227     }
1228 jgs 149
1229 jgs 82 bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1230     {
1231     return false;
1232     }
1233 jgs 104
1234 jgs 121 bool MeshAdapter::operator==(const AbstractDomain& other) const
1235 jgs 82 {
1236 jgs 121 const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
1237     if (temp!=0) {
1238     return (m_finleyMesh==temp->m_finleyMesh);
1239     } else {
1240     return false;
1241     }
1242 jgs 82 }
1243    
1244 jgs 121 bool MeshAdapter::operator!=(const AbstractDomain& other) const
1245 jgs 104 {
1246 jgs 121 return !(operator==(other));
1247 jgs 104 }
1248    
1249 jgs 150 int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const
1250 jgs 102 {
1251 jgs 150 int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1252     checkPasoError();
1253 jgs 102 return out;
1254     }
1255 jgs 149
1256 jgs 480 escript::Data MeshAdapter::getX() const
1257 jgs 102 {
1258     return continuousFunction(asAbstractContinuousDomain()).getX();
1259     }
1260 jgs 149
1261 jgs 480 escript::Data MeshAdapter::getNormal() const
1262 jgs 102 {
1263     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();
1264     }
1265 jgs 149
1266 jgs 480 escript::Data MeshAdapter::getSize() const
1267 jgs 102 {
1268     return function(asAbstractContinuousDomain()).getSize();
1269     }
1270    
1271 gross 1062 int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1272     {
1273     int *out=0,i;
1274     Finley_Mesh* mesh=m_finleyMesh.get();
1275     switch (functionSpaceType) {
1276     case(Nodes):
1277     if (mesh->Nodes!=NULL) {
1278     out=mesh->Nodes->Id;
1279     break;
1280     }
1281     case(ReducedNodes):
1282     throw FinleyAdapterException("Error - ReducedNodes not supported yet.");
1283     break;
1284     case(Elements):
1285     out=mesh->FaceElements->Id;
1286     break;
1287     case(ReducedElements):
1288     out=mesh->Elements->Id;
1289     break;
1290     case(FaceElements):
1291     out=mesh->FaceElements->Id;
1292     break;
1293     case(ReducedFaceElements):
1294     out=mesh->FaceElements->Id;
1295     break;
1296     case(Points):
1297     out=mesh->Points->Id;
1298     break;
1299     case(ContactElementsZero):
1300     out=mesh->ContactElements->Id;
1301     break;
1302     case(ReducedContactElementsZero):
1303     out=mesh->ContactElements->Id;
1304     break;
1305     case(ContactElementsOne):
1306     out=mesh->ContactElements->Id;
1307     break;
1308     case(ReducedContactElementsOne):
1309     out=mesh->ContactElements->Id;
1310     break;
1311     case(DegreesOfFreedom):
1312     out=mesh->Nodes->degreeOfFreedomId;
1313     break;
1314     case(ReducedDegreesOfFreedom):
1315     out=mesh->Nodes->reducedDegreeOfFreedomId;
1316     break;
1317     default:
1318     stringstream temp;
1319     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1320     throw FinleyAdapterException(temp.str());
1321     break;
1322     }
1323     return out;
1324     }
1325 jgs 110 int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1326     {
1327 gross 547 int out=0;
1328     Finley_Mesh* mesh=m_finleyMesh.get();
1329     switch (functionSpaceType) {
1330     case(Nodes):
1331     out=mesh->Nodes->Tag[sampleNo];
1332     break;
1333 gross 1062 case(ReducedNodes):
1334     throw FinleyAdapterException(" Error - ReducedNodes does not support tags.");
1335     break;
1336 gross 547 case(Elements):
1337     out=mesh->Elements->Tag[sampleNo];
1338     break;
1339 gross 1062 case(ReducedElements):
1340     out=mesh->Elements->Tag[sampleNo];
1341     break;
1342 gross 547 case(FaceElements):
1343     out=mesh->FaceElements->Tag[sampleNo];
1344     break;
1345 gross 1062 case(ReducedFaceElements):
1346     out=mesh->FaceElements->Tag[sampleNo];
1347     break;
1348 gross 547 case(Points):
1349     out=mesh->Points->Tag[sampleNo];
1350     break;
1351     case(ContactElementsZero):
1352     out=mesh->ContactElements->Tag[sampleNo];
1353     break;
1354 gross 1062 case(ReducedContactElementsZero):
1355     out=mesh->ContactElements->Tag[sampleNo];
1356     break;
1357 gross 547 case(ContactElementsOne):
1358     out=mesh->ContactElements->Tag[sampleNo];
1359     break;
1360 gross 1062 case(ReducedContactElementsOne):
1361     out=mesh->ContactElements->Tag[sampleNo];
1362 gross 547 break;
1363     case(DegreesOfFreedom):
1364 gross 1062 throw FinleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
1365 gross 547 break;
1366     case(ReducedDegreesOfFreedom):
1367 gross 1062 throw FinleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
1368 gross 547 break;
1369     default:
1370     stringstream temp;
1371     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1372     throw FinleyAdapterException(temp.str());
1373     break;
1374     }
1375     return out;
1376 jgs 110 }
1377    
1378 gross 1062
1379 gross 767 void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1380     {
1381     Finley_Mesh* mesh=m_finleyMesh.get();
1382     escriptDataC tmp=mask.getDataC();
1383     switch(functionSpaceType) {
1384     case(Nodes):
1385     Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
1386     break;
1387 gross 1062 case(ReducedNodes):
1388     throw FinleyAdapterException("Error - ReducedNodes does not support tags");
1389     break;
1390 gross 767 case(DegreesOfFreedom):
1391     throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
1392     break;
1393     case(ReducedDegreesOfFreedom):
1394     throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
1395     break;
1396     case(Elements):
1397     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1398     break;
1399 gross 1062 case(ReducedElements):
1400     Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
1401     break;
1402 gross 767 case(FaceElements):
1403     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1404     break;
1405 gross 1062 case(ReducedFaceElements):
1406     Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
1407     break;
1408 gross 767 case(Points):
1409     Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);
1410     break;
1411     case(ContactElementsZero):
1412     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1413     break;
1414 gross 1062 case(ReducedContactElementsZero):
1415     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1416     break;
1417 gross 767 case(ContactElementsOne):
1418     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1419     break;
1420 gross 1062 case(ReducedContactElementsOne):
1421     Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);
1422     break;
1423 gross 767 default:
1424     stringstream temp;
1425     temp << "Error - Finley does not know anything about function space type " << functionSpaceType;
1426     throw FinleyAdapterException(temp.str());
1427     }
1428     checkFinleyError();
1429     return;
1430     }
1431    
1432 gross 1044 void MeshAdapter::setTagMap(const std::string& name, int tag)
1433     {
1434     Finley_Mesh* mesh=m_finleyMesh.get();
1435     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
1436     checkPasoError();
1437     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
1438     }
1439 gross 767
1440 gross 1044 int MeshAdapter::getTag(const std::string& name) const
1441     {
1442     Finley_Mesh* mesh=m_finleyMesh.get();
1443     int tag=0;
1444     tag=Finley_Mesh_getTag(mesh, name.c_str());
1445     checkPasoError();
1446     // throwStandardException("MeshAdapter::getTag is not implemented.");
1447     return tag;
1448     }
1449    
1450     bool MeshAdapter::isValidTagName(const std::string& name) const
1451     {
1452     Finley_Mesh* mesh=m_finleyMesh.get();
1453     return Finley_Mesh_isValidTagName(mesh,name.c_str());
1454     }
1455    
1456     std::string MeshAdapter::showTagNames() const
1457     {
1458     stringstream temp;
1459     Finley_Mesh* mesh=m_finleyMesh.get();
1460     Finley_TagMap* tag_map=mesh->TagMap;
1461     while (tag_map) {
1462     temp << tag_map->name;
1463     tag_map=tag_map->next;
1464     if (tag_map) temp << ", ";
1465     }
1466     return temp.str();
1467     }
1468    
1469 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