/[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 751 - (hide annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years ago) by bcumming
File size: 39985 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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