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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26