/[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 682 - (hide annotations)
Mon Mar 27 02:43:09 2006 UTC (13 years, 6 months ago) by robwdcock
File size: 37983 byte(s)
+ NEW BUILD SYSTEM

This commit contains the new build system with cross-platform support.
Most things work are before though you can have more control.

ENVIRONMENT settings have changed:
+ You no longer require LD_LIBRARY_PATH or PYTHONPATH to point to the
esysroot for building and testing performed via scons
+ ACcESS altix users: It is recommended you change your modules to load
the latest intel compiler and other libraries required by boost to match
the setup in svn (you can override). The correct modules are as follows

module load intel_cc.9.0.026
export
MODULEPATH=${MODULEPATH}:/data/raid2/toolspp4/modulefiles/gcc-3.3.6
module load boost/1.33.0/python-2.4.1
module load python/2.4.1
module load numarray/1.3.3


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