/[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 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/finley/src/CPPAdapter/MeshAdapter.cpp
File size: 38304 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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