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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 817 - (hide annotations)
Sat Aug 26 03:08:52 2006 UTC (12 years, 7 months ago) by ksteube
Original Path: trunk/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 42220 byte(s)
Can now compile and run with MPI on shake71


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