Assembler.h 9.9 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
    virtual ~Assembler() 
    {}
73 74

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


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

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

    void calculateElementMatrix(const ElInfo *rowElInfo,
				const ElInfo *colElInfo,
91 92
				const ElInfo *smallElInfo,
				const ElInfo *largeElInfo,
93 94
				ElementMatrix *userMat,
				double factor = 1.0);
95 96 97 98

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

    /** \brief
     * Returns \ref rowFESpace.
     */
106 107
    inline const FiniteElemSpace* getRowFESpace() { 
      return rowFESpace; 
108
    }
109 110 111 112

    /** \brief
     * Returns \ref colFESpace.
     */
113 114
    inline const FiniteElemSpace* getColFESpace() { 
      return colFESpace; 
115
    }
116 117 118 119

    /** \brief
     * Returns \ref nRow.
     */
120 121
    inline int getNRow() { 
      return nRow; 
122
    }
123 124 125 126

    /** \brief
     * Returns \ref nCol.
     */
127 128
    inline int getNCol() { 
      return nCol; 
129
    }
130 131 132 133

    /** \brief
     * Sets \ref rememberElMat.
     */
134 135
    inline void rememberElementMatrix(bool rem) { 
      rememberElMat = rem; 
136
    }
137 138 139 140

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

    /** \brief
     * Returns \ref zeroOrderAssembler.
     */
    inline ZeroOrderAssembler* getZeroOrderAssembler() {
      return zeroOrderAssembler;
150
    }
151 152 153 154 155

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

    /** \brief
     * Returns \ref secondOrderAssembler.
     */
    inline SecondOrderAssembler* getSecondOrderAssembler() {
      return secondOrderAssembler;
168
    }
169 170 171 172

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

    /** \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) 
    {
192
      if (secondOrderAssembler) {
193 194 195 196
	TEST_EXIT(!secondOrderAssembler->getQuadrature())
	  ("quadrature already existing\n");
	secondOrderAssembler->setQuadrature(quad2);
      }
197
      if (firstOrderAssemblerGrdPsi) {
198 199 200 201
	TEST_EXIT(!firstOrderAssemblerGrdPsi->getQuadrature())
	  ("quadrature already existing\n");
	firstOrderAssemblerGrdPsi->setQuadrature(quad1GrdPsi);
      }
202
      if (firstOrderAssemblerGrdPhi) {
203 204 205 206
	TEST_EXIT(!firstOrderAssemblerGrdPhi->getQuadrature())
	  ("quadrature already existing\n");
	firstOrderAssemblerGrdPhi->setQuadrature(quad1GrdPhi);
      }
207
      if (zeroOrderAssembler) {
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
	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() {
225
      if (lastTraverseId != ElInfo::traverseId[omp_get_thread_num()]) {
226
	lastVecEl = lastMatEl = NULL;
227
	lastTraverseId = ElInfo::traverseId[omp_get_thread_num()];
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 268 269 270

    /** \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
271
    FirstOrderAssembler *firstOrderAssemblerGrdPsi;
272 273 274 275

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

    /** \brief
     * SubAssembler for the zero order terms
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    ZeroOrderAssembler *zeroOrderAssembler;
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 343 344 345

    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
346
    StandardAssembler(Operator *op,
347 348 349 350
		      Quadrature *quad2,
		      Quadrature *quad1GrdPsi,
		      Quadrature *quad1GrdPhi,
		      Quadrature *quad0,
Thomas Witkowski's avatar
Thomas Witkowski committed
351 352
		      const FiniteElemSpace *rowFESpace,
		      const FiniteElemSpace *colFESpace = NULL);
353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372
  };

  // ============================================================================
  // ===== 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
373
    OptimizedAssembler(Operator *op,
374 375 376 377
		       Quadrature *quad2,
		       Quadrature *quad1GrdPsi,
		       Quadrature *quad1GrdPhi,
		       Quadrature *quad0,
Thomas Witkowski's avatar
Thomas Witkowski committed
378 379
		       const FiniteElemSpace *rowFESpace,
		       const FiniteElemSpace *colFESpace = NULL);
380 381 382 383 384
  };

}

#endif // AMDIS_ASSEMBLER_H