SubAssembler.cc 6.47 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
101
102
103
104
105
				 Quadrature *quad)
  {
    // set corrdsAtQPs invalid
    coordsValid = false;

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

    // set gradients at QPs invalid
    std::map<const DOFVectorBase<double>*, GradientsAtQPs*>::iterator it2;
111
    for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2)
112
113
114
115
116
117
      ((*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) {
118
      if (largeElInfo == NULL)
Thomas Witkowski's avatar
Thomas Witkowski committed
119
	(*it)->initElement(smallElInfo, this, quad);
120
121
      else
	(*it)->initElement(smallElInfo, largeElInfo, this, quad);      
122
123
124
    }
  }

125

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

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

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

    // mark values as valid
    coordsValid = true;

    return coordsAtQPs;
  }

159

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

166
    TEST_EXIT_DBG(vec)("No DOF vector!\n");
167
168

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

    Quadrature *localQuad = quad ? quad : quadrature;

173
    if (!gradientsAtQPs[vec])
174
      gradientsAtQPs[vec] = new GradientsAtQPs;
175

176
177
    gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());

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

180
181
    const BasisFunction *psi = rowFeSpace->getBasisFcts();
    const BasisFunction *phi = colFeSpace->getBasisFcts();
182

183
    bool sameFeSpaces = 
184
      vec->getFeSpace() == rowFeSpace || vec->getFeSpace() == colFeSpace;
185

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

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

    return values;
  }

210

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

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

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

    Quadrature *localQuad = quad ? quad : quadrature;

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

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

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

236
}