17 |
|
|
18 |
#include <escript/FunctionSpace.h> |
#include <escript/FunctionSpace.h> |
19 |
#include <escript/FunctionSpaceFactory.h> |
#include <escript/FunctionSpaceFactory.h> |
20 |
|
|
21 |
|
#include <pasowrap/PatternBuilder.h> |
22 |
using namespace buckley; |
using namespace buckley; |
23 |
using namespace std; |
using namespace std; |
24 |
|
using namespace paso; |
25 |
|
|
26 |
|
|
27 |
namespace |
namespace |
415 |
throw BuckleyException("Not Implemented"); |
throw BuckleyException("Not Implemented"); |
416 |
} |
} |
417 |
|
|
418 |
|
|
419 |
|
|
420 |
|
|
421 |
|
|
422 |
escript::ASM_ptr BuckleyDomain::newSystemMatrix( |
escript::ASM_ptr BuckleyDomain::newSystemMatrix( |
423 |
const int row_blocksize, |
const int row_blocksize, |
424 |
const escript::FunctionSpace& row_functionspace, |
const escript::FunctionSpace& row_functionspace, |
426 |
const escript::FunctionSpace& column_functionspace, |
const escript::FunctionSpace& column_functionspace, |
427 |
const int type) const |
const int type) const |
428 |
{ |
{ |
429 |
throw BuckleyException("Not Implemented"); |
// This setup chunk copied from finley |
430 |
|
int reduceRowOrder=0; |
431 |
|
int reduceColOrder=0; |
432 |
|
// is the domain right? |
433 |
|
const BuckleyDomain& row_domain=dynamic_cast<const BuckleyDomain&>(*(row_functionspace.getDomain())); |
434 |
|
if (row_domain!=*this) |
435 |
|
throw BuckleyException("Error - domain of row function space does not match the domain of matrix generator."); |
436 |
|
const BuckleyDomain& col_domain=dynamic_cast<const BuckleyDomain&>(*(column_functionspace.getDomain())); |
437 |
|
if (col_domain!=*this) |
438 |
|
throw BuckleyException("Error - domain of columnn function space does not match the domain of matrix generator."); |
439 |
|
// is the function space type right |
440 |
|
|
441 |
|
// if (row_functionspace.getTypeCode()==DegreesOfFreedom) { |
442 |
|
reduceRowOrder=0; |
443 |
|
// } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) { |
444 |
|
// reduceRowOrder=1; |
445 |
|
// } else { |
446 |
|
// throw FinleyAdapterException("Error - illegal function space type for system matrix rows."); |
447 |
|
// } |
448 |
|
// if (column_functionspace.getTypeCode()==DegreesOfFreedom) { |
449 |
|
reduceColOrder=0; |
450 |
|
/* } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) { |
451 |
|
reduceColOrder=1;*/ |
452 |
|
// } else { |
453 |
|
// throw FinleyAdapterException("Error - illegal function space type for system matrix columns."); |
454 |
|
// } |
455 |
|
// generate matrix: |
456 |
|
|
457 |
|
// This is just here to remind me that I need to consider multiple ranks |
458 |
|
const unsigned int numranks=1; |
459 |
|
PatternBuilder* pb=makePB(numpts/numranks,26); |
460 |
|
|
461 |
|
// again, dummy values for a sole rank |
462 |
|
Paso_SystemMatrixPattern* psystemMatrixPattern=pb->generatePattern(0, numpts); |
463 |
|
|
464 |
|
|
465 |
|
// Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceRowOrder,reduceColOrder); |
466 |
|
// checkFinleyError(); |
467 |
|
// Paso_SystemMatrix* fsystemMatrix; |
468 |
|
// int trilinos = 0; |
469 |
|
// if (trilinos) { |
470 |
|
// #ifdef TRILINOS |
471 |
|
// /* Allocation Epetra_VrbMatrix here */ |
472 |
|
// #endif |
473 |
|
// } |
474 |
|
// else { |
475 |
|
// fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE); |
476 |
|
// } |
477 |
|
// checkPasoError(); |
478 |
|
// Paso_SystemMatrixPattern_free(fsystemMatrixPattern); |
479 |
|
// SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace); |
480 |
|
// return ASM_ptr(sma); |
481 |
|
|
482 |
|
|
483 |
|
|
484 |
|
throw BuckleyException("Not Implemented ::newSystemMatrix"); |
485 |
} |
} |
486 |
|
|
487 |
escript::ATP_ptr BuckleyDomain::newTransportProblem( |
escript::ATP_ptr BuckleyDomain::newTransportProblem( |
748 |
if (j>3){lz*=3;} |
if (j>3){lz*=3;} |
749 |
if (j==1 || j==2 || j==5 || j==6) {lx*=3;} |
if (j==1 || j==2 || j==5 || j==6) {lx*=3;} |
750 |
if (j==2 || j==3 || j==6 || j==7) {ly*=3;} |
if (j==2 || j==3 || j==6 || j==7) {ly*=3;} |
751 |
|
|
752 |
|
// is an element with a kink in it a problem here? |
753 |
|
|
754 |
dest[getRelIndex(MAT3X3, 0, 0)]=values2[j*3+0]/lx; |
dest[getRelIndex(MAT3X3, 0, 0)]=values2[j*3+0]/lx; |
755 |
dest[getRelIndex(MAT3X3, 0, 1)]=values2[j*3+0]/ly; |
dest[getRelIndex(MAT3X3, 0, 1)]=values2[j*3+0]/ly; |
756 |
dest[getRelIndex(MAT3X3, 0, 2)]=values2[j*3+0]/lz; |
dest[getRelIndex(MAT3X3, 0, 2)]=values2[j*3+0]/lz; |