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