Assembler.h 9.85 KB
Newer Older
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
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file Assembler.h */

/**
 * \defgroup Assembler Assembler module
 *
 * \brief
 * Contains the operator and assembler classes:
 * @{ <img src="assembler.png"> @}
 */

#ifndef AMDIS_ASSEMBLER_H
#define AMDIS_ASSEMBLER_H

#include <vector>
#include "FixVec.h"
#include "MemoryManager.h"
36 37 38 39
#include "ZeroOrderAssembler.h"
#include "FirstOrderAssembler.h"
#include "SecondOrderAssembler.h"
#include "ElInfo.h"
40
#include "OpenMP.h"
41 42 43

namespace AMDiS {

44
  //  class ElInfo;
45 46 47 48
  class Element;
  class Quadrature;
  class ElementMatrix;
  class ElementVector;
49
  class Operator;
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

  /**
   * \ingroup Assembler
   * 
   * \brief
   * Assembles element matrices and vectors for a given Operator. Uses
   * one SubAssembler for all second order terms of the Operator, one for all
   * first order terms, and one for all zero order terms.
   */
  class Assembler
  {
  public:
    MEMORY_MANAGED(Assembler);

    /** \brief
     * Constructor.
     */
67 68 69
    Assembler(Operator *op,
	      const FiniteElemSpace *rowFESpace,
	      const FiniteElemSpace *colFESpace = NULL);
70 71 72 73

    virtual ~Assembler() {};

    ElementMatrix *initElementMatrix(ElementMatrix *elMat, 
74 75
				     const ElInfo *rowElInfo,
				     const ElInfo *colElInfo = NULL);
76 77 78 79 80 81 82 83


    ElementVector *initElementVector(ElementVector *elVec,
				     const ElInfo *elInfo);

    /** \brief
     * Assembles the element matrix for the given ElInfo
     */
84 85 86 87 88 89 90 91
    void calculateElementMatrix(const ElInfo *elInfo, 
				ElementMatrix *userMat, 
				double factor = 1.0);

    void calculateElementMatrix(const ElInfo *rowElInfo,
				const ElInfo *colElInfo,
				ElementMatrix *userMat,
				double factor = 1.0);
92 93 94 95

    /** \brief
     * Assembles the element vector for the given ElInfo
     */
96 97 98
    void calculateElementVector(const ElInfo *elInfo, 
				ElementVector *userVec, 
				double factor = 1.0);
99 100 101 102

    /** \brief
     * Returns \ref rowFESpace.
     */
103 104 105
    inline const FiniteElemSpace* getRowFESpace() { 
      return rowFESpace; 
    };
106 107 108 109

    /** \brief
     * Returns \ref colFESpace.
     */
110 111 112
    inline const FiniteElemSpace* getColFESpace() { 
      return colFESpace; 
    };
113 114 115 116

    /** \brief
     * Returns \ref nRow.
     */
117 118 119
    inline int getNRow() { 
      return nRow; 
    };
120 121 122 123

    /** \brief
     * Returns \ref nCol.
     */
124 125 126
    inline int getNCol() { 
      return nCol; 
    };
127 128 129 130

    /** \brief
     * Sets \ref rememberElMat.
     */
131 132 133
    inline void rememberElementMatrix(bool rem) { 
      rememberElMat = rem; 
    };
134 135 136 137

    /** \brief
     * Sets \ref rememberElVec.
     */
138 139 140
    inline void rememberElementVector(bool rem) { 
      rememberElVec = rem; 
    };
141 142 143 144 145 146 147 148 149 150 151 152

    /** \brief
     * Returns \ref zeroOrderAssembler.
     */
    inline ZeroOrderAssembler* getZeroOrderAssembler() {
      return zeroOrderAssembler;
    };

    /** \brief
     * Returns \ref firstOrderAssemblerGrdPsi or \ref firstOrderAssemblerGrdPhi
     * depending on type.
     */
153
    inline FirstOrderAssembler* getFirstOrderAssembler(FirstOrderType type = GRD_PSI) 
154
    {
155
      return (type == GRD_PSI) ? 
156 157 158 159 160 161 162 163 164 165 166 167 168 169
	firstOrderAssemblerGrdPsi : 
	firstOrderAssemblerGrdPhi;
    };

    /** \brief
     * Returns \ref secondOrderAssembler.
     */
    inline SecondOrderAssembler* getSecondOrderAssembler() {
      return secondOrderAssembler;
    };

    /** \brief
     * Returns \ref operat;
     */
170 171 172
    inline Operator* getOperator() { 
      return operat; 
    };
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188

    /** \brief
     * Initialisation for the given ElInfo. The call is deligated to
     * the sub assemblers.
     */
    void initElement(const ElInfo *elInfo,
		     Quadrature *quad = NULL);

    /** \brief
     * Sets quadratures of all sub assemblers.
     */
    void setQuadratures(Quadrature *quad2,
			Quadrature *quad1GrdPsi,
			Quadrature *quad1GrdPhi,
			Quadrature *quad0) 
    {
189
      if (secondOrderAssembler) {
190 191 192 193
	TEST_EXIT(!secondOrderAssembler->getQuadrature())
	  ("quadrature already existing\n");
	secondOrderAssembler->setQuadrature(quad2);
      }
194
      if (firstOrderAssemblerGrdPsi) {
195 196 197 198
	TEST_EXIT(!firstOrderAssemblerGrdPsi->getQuadrature())
	  ("quadrature already existing\n");
	firstOrderAssemblerGrdPsi->setQuadrature(quad1GrdPsi);
      }
199
      if (firstOrderAssemblerGrdPhi) {
200 201 202 203
	TEST_EXIT(!firstOrderAssemblerGrdPhi->getQuadrature())
	  ("quadrature already existing\n");
	firstOrderAssemblerGrdPhi->setQuadrature(quad1GrdPhi);
      }
204
      if (zeroOrderAssembler) {
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
	TEST_EXIT(!zeroOrderAssembler->getQuadrature())
	  ("quadrature already existing\n");
	zeroOrderAssembler->setQuadrature(quad0);
      }
    };

  protected:
    /** \brief
     * Vector assembling by element matrix-vector multiplication. 
     * Usefull if an element matrix was already calculated.
     */
    void matVecAssemble(const ElInfo *elInfo, ElementVector *vec);

    /** \brief
     * Checks whether this is a new travese.
     */
    inline void checkForNewTraverse() {
222
      if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) {
223
	lastVecEl = lastMatEl = NULL;
224
	lastTraverseId = ElInfo::traverseId[omp_get_thread_num()];
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
      }
    };

    /** \brief
     * Checks whether quadratures for sub assemblers are already set.
     * If not they will be created.
     */
    void checkQuadratures();

  protected:
    /** \brief
     * Operator this Assembler belongs to.
     */
    Operator *operat;

    /** \brief
     * Row FiniteElemSpace.
     */
    const FiniteElemSpace *rowFESpace;

    /** \brief
     * Column FiniteElemSpace.
     */
    const FiniteElemSpace *colFESpace;

    /** \brief
     * Number of rows.
     */
    int nRow;

    /** \brief
     * Number of columns.
     */
    int nCol;

    /** \brief
     * SubAssembler for the second order terms
     */
    SecondOrderAssembler *secondOrderAssembler;

    /** \brief
     * SubAssembler for the first order terms (grdPsi)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
268
    FirstOrderAssembler *firstOrderAssemblerGrdPsi;
269 270 271 272

    /** \brief
     * SubAssembler for the first order terms (grdPhi)
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    FirstOrderAssembler *firstOrderAssemblerGrdPhi;  
274 275 276 277

    /** \brief
     * SubAssembler for the zero order terms
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
278
    ZeroOrderAssembler *zeroOrderAssembler;
279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342

    bool remember;

    /** \brief
     * Determines whether the element matrix should be stored locally.
     */
    bool rememberElMat;

    /** \brief
     * Determines whether the element vector should be stored locally.
     */
    bool rememberElVec;

    /** \brief
     * locally stored element matrix
     */
    ElementMatrix* elementMatrix;

    /** \brief
     * locally stored element vector
     */
    ElementVector* elementVector;
  
    /** \brief
     * Used to check whether \ref initElement() must be called, because
     * a new Element is visited.
     */
    Element* lastMatEl;

    /** \brief
     * Used to check whether \ref initElement() must be called, because
     * a new Element is visited.
     */
    Element* lastVecEl;

    /** \brief
     * Used to check for new traverse.
     */
    int lastTraverseId;

    friend class SubAssembler;
    friend class ZeroOrderAssembler;
    friend class FirstOrderAssembler;
    friend class SecondOrderAssembler;
  };

  // ============================================================================
  // ===== class StandardAssembler =============================================
  // ============================================================================

  /**
   * \ingroup Assembler
   *
   * \brief
   * Assembler using non optimized sub assemblers.
   */
  class StandardAssembler : public Assembler
  {
  public:
    MEMORY_MANAGED(StandardAssembler);

    /** \brief
     * Constructor.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    StandardAssembler(Operator *op,
344 345 346 347
		      Quadrature *quad2,
		      Quadrature *quad1GrdPsi,
		      Quadrature *quad1GrdPhi,
		      Quadrature *quad0,
Thomas Witkowski's avatar
Thomas Witkowski committed
348 349
		      const FiniteElemSpace *rowFESpace,
		      const FiniteElemSpace *colFESpace = NULL);
350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369
  };

  // ============================================================================
  // ===== class OptimizedAssembler =============================================
  // ============================================================================

  /**
   * \ingroup Assembler
   *
   * \brief
   * Assembler using optimized sub assemblers.
   */
  class OptimizedAssembler : public Assembler
  {
  public:
    MEMORY_MANAGED(OptimizedAssembler);

    /** \brief
     * Constructor.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
370
    OptimizedAssembler(Operator *op,
371 372 373 374
		       Quadrature *quad2,
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0,
Thomas Witkowski's avatar
Thomas Witkowski committed
375 376
		       const FiniteElemSpace *rowFESpace,
		       const FiniteElemSpace *colFESpace = NULL);
377 378 379 380 381
  };

}

#endif // AMDIS_ASSEMBLER_H