/[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 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (13 years, 9 months ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 37731 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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