/[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 798 - (hide annotations)
Fri Aug 4 01:05:36 2006 UTC (13 years, 3 months ago) by gross
File size: 42180 byte(s)
Reimplementation of the assemblage with persistent jacobeans.
There are also a few changes to the tests which has now
dramatically reduced the memory demand.


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