SubAssembler.cc 6.52 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include "SubAssembler.h"
#include "Assembler.h"
#include "FiniteElemSpace.h"
#include "Operator.h"
#include "BasisFunction.h"
#include "OpenMP.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "DOFVector.h"

namespace AMDiS {

  SubAssembler::SubAssembler(Operator *op,
			     Assembler *assembler,
			     Quadrature *quadrat,
			     int order, 
			     bool optimized,
			     FirstOrderType type) 
19
20
21
    : rowFeSpace(assembler->rowFeSpace),
      colFeSpace(assembler->colFeSpace),
      nRow(0),
22
23
24
25
26
27
28
29
      nCol(0),
      coordsAtQPs(NULL),
      coordsNumAllocated(0),
      quadrature(quadrat),
      psiFast(NULL),
      phiFast(NULL),
      symmetric(true),
      opt(optimized),
30
31
      firstCall(true),
      name("")
32
  {
33
34
35
36
37
38
39
    FUNCNAME("SubAssembler::SubAssembler()");

    TEST_EXIT(rowFeSpace)("No row FE space defined!\n");
    TEST_EXIT(colFeSpace)("No col FE space defined!\n");

    const BasisFunction *psi = rowFeSpace->getBasisFcts();
    const BasisFunction *phi = colFeSpace->getBasisFcts();
40
41
42
43

    nRow = psi->getNumber();
    nCol = phi->getNumber();

44
    int maxThreads = omp_get_overall_max_threads();
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
    terms.resize(maxThreads);

    switch (order) {
    case 0:
      terms = op->zeroOrder;
      break;
    case 1:
      if (type == GRD_PHI)
	terms = op->firstOrderGrdPhi;
      else 
	terms = op->firstOrderGrdPsi;
      break;
    case 2:
      terms = op->secondOrder;
      break;
    }

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
    // If the number of basis functions in row and col fe space are equal,
    // the element matrix may be symmetric if all operator terms are symmetric.
    // If the row and col fe space are different, the element matrix is neither
    // quadratic, and therefore cannot be symmetric.
    if (nRow == nCol) {
      symmetric = true;
      for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
	if (!(terms[0][i])->isSymmetric()) {
	  symmetric = false;
	  break;
	}
      }  
    } else {
      symmetric = false;
    }
77

78
    dim = rowFeSpace->getMesh()->getDim();
79
80
  }

81

82
83
84
85
86
  FastQuadrature *SubAssembler::updateFastQuadrature(FastQuadrature *quadFast,
						     const BasisFunction *psi,
						     Flag updateFlag)
  {
    if (!quadFast) {
Thomas Witkowski's avatar
Thomas Witkowski committed
87
      quadFast = FastQuadrature::provideFastQuadrature(psi, *quadrature, updateFlag);
88
89
90
91
92
93
94
95
    } else {
      if (!quadFast->initialized(updateFlag))
	quadFast->init(updateFlag);
    }

    return quadFast;
  }

96

Thomas Witkowski's avatar
Thomas Witkowski committed
97
98
  void SubAssembler::initElement(const ElInfo* smallElInfo,
				 const ElInfo *largeElInfo,
99
100
				 Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
    FUNCNAME("SubAssembler::initElement()");

103
104
105
106
107
    // set corrdsAtQPs invalid
    coordsValid = false;

    // set values at QPs invalid
    std::map<const DOFVectorBase<double>*, ValuesAtQPs*>::iterator it1;
108
    for (it1 = valuesAtQPs.begin(); it1 != valuesAtQPs.end(); ++it1)
109
110
111
112
      ((*it1).second)->valid = false;

    // set gradients at QPs invalid
    std::map<const DOFVectorBase<double>*, GradientsAtQPs*>::iterator it2;
113
    for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2)
114
115
116
117
118
119
      ((*it2).second)->valid = false;

    int myRank = omp_get_thread_num();
    // calls initElement of each term
    std::vector<OperatorTerm*>::iterator it;
    for (it = terms[myRank].begin(); it != terms[myRank].end(); ++it) {
120
      if (largeElInfo == NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
121
	(*it)->initElement(smallElInfo, this, quad);
122
123
      else
	(*it)->initElement(smallElInfo, largeElInfo, this, quad);      
124
125
126
    }
  }

127

Thomas Witkowski's avatar
Thomas Witkowski committed
128
  WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo, 
129
130
131
132
133
134
135
						    Quadrature *quad)
  {
    Quadrature *localQuad = quad ? quad : quadrature;
  
    const int nPoints = localQuad->getNumPoints();

    // already calculated for this element ?
Thomas Witkowski's avatar
Thomas Witkowski committed
136
    if (coordsValid)
137
138
139
140
      return coordsAtQPs;
   
    if (coordsAtQPs)  {
      if (coordsNumAllocated != nPoints) {
Thomas Witkowski's avatar
Thomas Witkowski committed
141
142
	delete [] coordsAtQPs;
        coordsAtQPs = new WorldVector<double>[nPoints];
143
144
145
	coordsNumAllocated = nPoints;
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
146
      coordsAtQPs = new WorldVector<double>[nPoints];
147
148
149
150
151
      coordsNumAllocated = nPoints;
    }

    // set new values
    WorldVector<double>* k = &(coordsAtQPs[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
152
    for (int l = 0; k < &(coordsAtQPs[nPoints]); ++k, ++l)
153
      elInfo->coordToWorld(localQuad->getLambda(l), *k);
154
155
156
157
158
159
160

    // mark values as valid
    coordsValid = true;

    return coordsAtQPs;
  }

161

162
  WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* vec, 
163
164
165
166
167
						       const ElInfo* elInfo,
						       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getGradientsAtQPs()");

168
    TEST_EXIT_DBG(vec)("No DOF vector!\n");
169
170

    if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid) 
Thomas Witkowski's avatar
Thomas Witkowski committed
171
      return &(gradientsAtQPs[vec]->values[0]);
172
173
174

    Quadrature *localQuad = quad ? quad : quadrature;

175
    if (!gradientsAtQPs[vec])
176
      gradientsAtQPs[vec] = new GradientsAtQPs;
177

178
179
    gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());

Thomas Witkowski's avatar
Thomas Witkowski committed
180
    WorldVector<double> *values = &(gradientsAtQPs[vec]->values[0]);
181

182
183
    const BasisFunction *psi = rowFeSpace->getBasisFcts();
    const BasisFunction *phi = colFeSpace->getBasisFcts();
184

185
    bool sameFeSpaces = 
186
      vec->getFeSpace() == rowFeSpace || vec->getFeSpace() == colFeSpace;
187

188
189
    if (opt && !quad && sameFeSpaces) {
      if (vec->getFeSpace()->getBasisFcts() == psi)
190
	psiFast = updateFastQuadrature(psiFast, psi, INIT_GRD_PHI);
191
      else if(vec->getFeSpace()->getBasisFcts() == phi)
192
193
194
195
	phiFast = updateFastQuadrature(phiFast, phi, INIT_GRD_PHI);
    }
  
    // calculate new values
196
    const BasisFunction *basFcts = vec->getFeSpace()->getBasisFcts();
197

198
    if (opt && !quad && sameFeSpaces) {
199
      if (psiFast->getBasisFunctions() == basFcts)
200
	vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
201
202
      else
	vec->getGrdAtQPs(elInfo, NULL, phiFast, values);      
203
204
205
206
207
208
209
210
211
    } else {
      vec->getGrdAtQPs(elInfo, localQuad, NULL, values);
    }
   
    gradientsAtQPs[vec]->valid = true;

    return values;
  }

212

213
  WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
217
218
219
						       const ElInfo* smallElInfo,
						       const ElInfo* largeElInfo,
						       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getGradientsAtQPs()");

220
    TEST_EXIT_DBG(vec)("No DOF vector!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
221
222

    if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid) 
Thomas Witkowski's avatar
Thomas Witkowski committed
223
      return &(gradientsAtQPs[vec]->values[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
224
225
226

    Quadrature *localQuad = quad ? quad : quadrature;

227
228
    if (!gradientsAtQPs[vec])
      gradientsAtQPs[vec] = new GradientsAtQPs;    
Thomas Witkowski's avatar
Thomas Witkowski committed
229
230
    gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());

Thomas Witkowski's avatar
Thomas Witkowski committed
231
    WorldVector<double> *values = &(gradientsAtQPs[vec]->values[0]);
Thomas Witkowski's avatar
Thomas Witkowski committed
232
233
234
235
236
237

    gradientsAtQPs[vec]->valid = true;
    vec->getGrdAtQPs(smallElInfo, largeElInfo, localQuad, NULL, values);
    return values;
  }

238
}