parametric.cc 5.04 KB
 Peter Gottschling committed Feb 15, 2008 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 #include "AMDiS.h" using namespace std; using namespace AMDiS; // =========================================================================== // ===== function definitions ================================================ // =========================================================================== /** \brief * RHS function */ class F : public AbstractFunction > { public: MEMORY_MANAGED(F); F(int degree) : AbstractFunction >(degree) {}; /** \brief * Implementation of AbstractFunction::operator(). */ const double& operator()(const WorldVector& x) const { static double result = 0.0; result = -2 * x[0]; return result; }; }; /** \brief * Implementation of a scalar problem. In \ref buildBeforeCoarsen() parametric * coordinates for the vertices are filled in a DOFVector. This DOFVector is * used in \ref parametric to parametrize elements while mesh traversal. */ class ParametricSphere : public ProblemScal { public: /** \brief * constructor */ ParametricSphere(const char *name) : ProblemScal(name), parametric(NULL) {}; /** \brief * destructor */ ~ParametricSphere() { int i; for(i = 0; i < Global::getGeo(WORLD); i++) { DELETE parametricCoords[i]; } DELETE parametric; } /** \brief * initialization of the base class and creation of \ref parametric. */ void initialize(Flag initFlag, ProblemScal *adoptProblem = NULL, Flag adoptFlag = INIT_NOTHING) { ProblemScal::initialize(initFlag, adoptProblem, adoptFlag); // ===== create projection ===== WorldVector ballCenter; ballCenter.set(0.0); NEW BallProject(1, VOLUME_PROJECTION, ballCenter, 1.0); // ===== create coords vector ===== int i; for(i = 0; i < Global::getGeo(WORLD); i++) { parametricCoords[i] = NEW DOFVector(this->getFESpace(), "parametric coords"); } // ===== create parametric object ===== parametric = NEW ParametricFirstOrder(¶metricCoords); // ===== enable parametric traverse ===== this->getMesh()->setParametric(parametric); }; /** \brief * Implementation of ProblemStatBase::buildBeforeCoarsen(). */ void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag flag) { FUNCNAME("ParametricSphere::buildAfterCoarsen()"); MSG("calculation of parametric coordinates\n"); int i, j; int preDOFs = this->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX); int dim = this->getMesh()->getDim(); int dow = Global::getGeo(WORLD); WorldVector vertexCoords; const DegreeOfFreedom **dof; DegreeOfFreedom vertexIndex; // ===== disable parametric traverse ===== this->getMesh()->setParametric(NULL); TraverseStack stack; ElInfo *elInfo = NULL; elInfo = stack.traverseFirst(this->getMesh(), -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS); while(elInfo) { dof = elInfo->getElement()->getDOF(); for(i = 0; i < dim + 1; i++) { vertexCoords = elInfo->getCoord(i); Projection::getProjection(1)->project(vertexCoords); for(j = 0; j < dow; j++) { (*(parametricCoords[j]))[dof[i][preDOFs]] = vertexCoords[j]; } } elInfo = stack.traverseNext(elInfo); } // ===== enable parametric traverse ===== this->getMesh()->setParametric(parametric); }; protected: /** \brief * DOFVector storing parametric coordinates. */ WorldVector*> parametricCoords; /** \brief * Parametrizes vertex coordinates while mesh traversal. */ ParametricFirstOrder *parametric; }; // =========================================================================== // ===== main program ======================================================== // =========================================================================== int main(int argc, char* argv[]) { FUNCNAME("parametric main"); // ===== check for init file ===== TEST_EXIT(argc == 2)("usage: parametric initfile\n"); // ===== init parameters ===== Parameters::init(false, argv[1]); // ===== create and init the scalar problem ===== ParametricSphere parametric("parametric"); parametric.initialize(INIT_ALL); // === create adapt info === AdaptInfo *adaptInfo = NEW AdaptInfo("parametric->adapt", 1); // === create adapt === AdaptStationary *adapt = NEW AdaptStationary("parametric->adapt", ¶metric, adaptInfo); // ===== create matrix operator ===== Operator matrixOperator(Operator::MATRIX_OPERATOR, parametric.getFESpace()); matrixOperator.addSecondOrderTerm(NEW Laplace_SOT); parametric.addMatrixOperator(&matrixOperator); // ===== create rhs operator ===== Operator rhsOperator(Operator::VECTOR_OPERATOR, parametric.getFESpace()); int degree = parametric.getFESpace()->getBasisFcts()->getDegree(); rhsOperator.addZeroOrderTerm(NEW CoordsAtQP_ZOT(NEW F(degree))); parametric.addVectorOperator(&rhsOperator); // ===== start adaption loop ===== adapt->adapt(); parametric.writeFiles(adaptInfo, true); return 0; }