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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 624 - (hide annotations)
Wed Mar 22 23:00:01 2006 UTC (16 years, 10 months ago) by robwdcock
File size: 37989 byte(s) 
+ Finish off the header file path changes
+ The C/C++ libraries now all compile on the ivec altix (cognac) using the xplatform build system
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 624 #include "escriptcpp/Data.h"
19     #include "escriptcpp/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