Assembler.cc 45.8 KB
Newer Older
1001

1002
    for (int iq = 0; iq < numPoints; iq++) {
1003
1004
      Lb[iq].set(0.0);
    }
1005
1006
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
1007
1008
    }
  
1009
    for (int iq = 0; iq < numPoints; iq++) {
1010
1011
1012
1013

      Lb[iq] *= elInfo->getDet();
      grdPsi = psiFast->getGradient(iq);

1014
      for (int i = 0; i < nRow; i++) {
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
	(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
      }
    }
  }

  Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
  {
  }

  void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);

    const int *l;
    const double *values;

1032
    if (firstCall) {
1033
1034
1035
1036
1037
1038
1039
1040
1041
      q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int **nEntries = q01->getNumberEntries();
1042
    int myRank = omp_get_thread_num();
1043
    Lb[0].set(0.0);
1044
1045
1046

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
1047
1048
1049
1050
    }

    Lb[0] *= elInfo->getDet();

1051
1052
1053
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	l = q01->getLVec(i, j);
1054
	values = q01->getValVec(i, j);
1055
	double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1056
1057
1058
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][l[m]];
	}
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
	(*mat)[i][j] += val;
      }
    }
  }

  void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    VectorOfFixVecs<DimVec<double> > Lb(dim,1,NO_INIT);

    const int *k;
    const double *values;

Thomas Witkowski's avatar
Thomas Witkowski committed
1071
    if (firstCall) {
1072
1073
1074
1075
1076
1077
1078
1079
1080
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int *nEntries = q1->getNumberEntries();
1081
    int myRank = omp_get_thread_num();
1082
    Lb[0].set(0.0);
1083
1084
1085

    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>(terms[myRank][i]))->getLb(elInfo, 1, Lb);
1086
1087
1088
1089
    }

    Lb[0] *= elInfo->getDet();

Thomas Witkowski's avatar
Thomas Witkowski committed
1090
1091
    for (int i = 0; i < nRow; i++) {
      k = q1->getKVec(i);
1092
      values = q1->getValVec(i);
1093
      double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1094
1095
1096
      for (int m = 0; m < nEntries[i]; m++) {
	val += values[m] * Lb[0][k[m]];
      }
1097
1098
1099
1100
1101
1102
      (*vec)[i] += val;
    }
  }

  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
  {
    q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
				      owner->getColFESpace()->getBasisFcts(), 
				      quadrature);
    tmpLALt.resize(omp_get_max_threads());
    for (int i = 0; i < omp_get_max_threads(); i++) {
      tmpLALt[i] = NEW DimMat<double>*;
      *(tmpLALt[i]) = NEW DimMat<double>(dim, NO_INIT);
    }
  }

  Pre2::~Pre2()
  {
    for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
      DELETE *(tmpLALt[i]);
      DELETE tmpLALt[i];
    }
  }
1121
1122
1123
1124
1125
1126
1127

  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    const int **nEntries;
    const int *k, *l;
    const double *values;

1128
    int myRank = omp_get_thread_num();
1129
1130
1131
    DimMat<double> **LALt = tmpLALt[myRank];
    DimMat<double> &tmpMat = *LALt[0];
    tmpMat.set(0.0);
1132

1133
1134
    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
1135
1136
    }

1137
    tmpMat *= elInfo->getDet();
1138
1139
1140
1141

    nEntries = q11->getNumberEntries();

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1142
1143
1144
      for (int i = 0; i < nRow; i++) {
	k = q11->getKVec(i, i);
	l = q11->getLVec(i, i);
1145
	values = q11->getValVec(i, i);
1146
	double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1147
	for (int m = 0; m < nEntries[i][i]; m++) {
1148
	  val += values[m] * tmpMat[k[m]][l[m]];
Thomas Witkowski's avatar
Thomas Witkowski committed
1149
	}
1150
	(*mat)[i][i] += val;
Thomas Witkowski's avatar
Thomas Witkowski committed
1151
1152
1153
1154

	for (int j = i + 1; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
1155
	  values = q11->getValVec(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
1156
1157
	  val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
1158
	    val += values[m] * tmpMat[k[m]][l[m]];
Thomas Witkowski's avatar
Thomas Witkowski committed
1159
	  }
1160
1161
1162
1163
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1164
1165
1166
1167
1168
    } else {  /*  A not symmetric or psi != phi        */
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
1169
	  values = q11->getValVec(i, j);
1170
	  double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1171
	  for (int m = 0; m < nEntries[i][j]; m++) {
1172
	    val += values[m] * tmpMat[k[m]][l[m]];
Thomas Witkowski's avatar
Thomas Witkowski committed
1173
	  }
1174
1175
1176
1177
1178
1179
1180
1181
	  (*mat)[i][j] += val;
	}
      }
    }
  }

  Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
Thomas Witkowski's avatar
Thomas Witkowski committed
1182
  {}
1183
1184
1185
1186
1187
1188

  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    double val;
    VectorOfFixVecs<DimVec<double> > *grdPsi, *grdPhi;

Thomas Witkowski's avatar
Thomas Witkowski committed
1189
    if (firstCall) {
1190
1191
1192
1193
1194
1195
1196
1197
1198
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
      firstCall = false;
    }

    int nPoints = quadrature->getNumPoints();
    DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
1199
1200
    int myRank = omp_get_thread_num();

Thomas Witkowski's avatar
Thomas Witkowski committed
1201
1202
1203
1204
    for (int i = 0; i < nPoints;i++) {
      LALt[i] = NEW DimMat<double>(dim, NO_INIT);
    }
    for (int iq = 0; iq < nPoints; iq++) {
1205
1206
      LALt[iq]->set(0.0);
    }
1207
1208
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
1209
1210
1211
    }

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1212
      for (int iq = 0; iq < nPoints; iq++) {
1213
1214
1215
1216
1217
	(*LALt[iq]) *= elInfo->getDet();

	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

Thomas Witkowski's avatar
Thomas Witkowski committed
1218
	for (int i = 0; i < nRow; i++) {
1219
1220
1221
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
	    ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));

Thomas Witkowski's avatar
Thomas Witkowski committed
1222
	  for (int j = i + 1; j < nCol; j++) {
1223
1224
1225
1226
1227
1228
1229
1230
	    val = quadrature->getWeight(iq) * ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    }
    else {      /*  non symmetric assembling   */
Thomas Witkowski's avatar
Thomas Witkowski committed
1231
      for (int iq = 0; iq < nPoints; iq++) {
1232
1233
1234
1235
1236
	(*LALt[iq]) *= elInfo->getDet();

	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

Thomas Witkowski's avatar
Thomas Witkowski committed
1237
1238
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
1239
1240
1241
1242
1243
1244
1245
	    (*mat)[i][j] += quadrature->getWeight(iq) *
	      ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
	  }
	}
      }
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1246
1247
1248
    for (int i = 0;i < nPoints; i++) 
      DELETE LALt[i];

1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
    DELETE [] LALt; 
  }

  Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, false)
  {}

  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    double val;
    DimVec<double> grdPsi(dim, NO_INIT);
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();

    DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
Thomas Witkowski's avatar
Thomas Witkowski committed
1268
1269
    for (int iq = 0; iq < nPoints; iq++) {
      LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
1270
1271
      LALt[iq]->set(0.0);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1272

1273
1274
1275
1276
    int myRank = omp_get_thread_num();

    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
1277
1278
1279
    }

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1280
      for (int iq = 0; iq < nPoints; iq++) {
1281
1282
	(*LALt[iq]) *= elInfo->getDet();

Thomas Witkowski's avatar
Thomas Witkowski committed
1283
	for (int i = 0; i < nCol; i++) {
1284
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
1285
1286
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
1287
	for (int i = 0; i < nRow; i++) {
1288
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
1289
1290
1291
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));

Thomas Witkowski's avatar
Thomas Witkowski committed
1292
	  for (int j = i + 1; j < nCol; j++) {
1293
1294
1295
1296
1297
1298
	    val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1299
1300
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
1301
1302
	(*LALt[iq]) *= elInfo->getDet();

Thomas Witkowski's avatar
Thomas Witkowski committed
1303
	for (int i = 0; i < nCol; i++) {
1304
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
1305
1306
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
1307
	for (int i = 0; i < nRow; i++) {
1308
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
Thomas Witkowski's avatar
Thomas Witkowski committed
1309
	  for (int j = 0; j < nCol; j++) {
1310
1311
1312
1313
1314
1315
1316
	    (*mat)[i][j] += quadrature->getWeight(iq) *
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	  }
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1317
1318
1319
    for (int iq = 0; iq < nPoints; iq++) 
      DELETE LALt[iq];

1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
    DELETE [] LALt;
  }

  Assembler::Assembler(Operator  *op,
		       const FiniteElemSpace *rowFESpace_,
		       const FiniteElemSpace *colFESpace_) 
    : operat(op),
      rowFESpace(rowFESpace_),
      colFESpace(colFESpace_ ? colFESpace_ : rowFESpace_),
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
      elementMatrix(NULL),
      elementVector(NULL),
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1340
  {}
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

    Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
1357
    
1358
1359
1360
    checkForNewTraverse();

    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
1361
    
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(elInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
  
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

    if (secondOrderAssembler)
      secondOrderAssembler->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPhi)
      firstOrderAssemblerGrdPhi->calculateElementMatrix(elInfo, mat);
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementMatrix(elInfo, mat);

Thomas Witkowski's avatar
Thomas Witkowski committed
1389
    if (rememberElMat && userMat) {
1390
      axpy(factor, *elementMatrix, *userMat);
1391
    }      
1392
1393
1394
1395
1396
1397
1398
1399
  }

  void Assembler::calculateElementVector(const ElInfo *elInfo, 
					 ElementVector *userVec,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
1400
    if (remember && factor != 1.0) {
1401
1402
1403
      rememberElVec = true;
    }

1404
    if (rememberElVec && !elementVector)
1405
1406
1407
1408
1409
1410
1411
1412
      elementVector = NEW ElementVector(nRow);

    Element *el = elInfo->getElement();

    checkForNewTraverse();

    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
1413
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
1414
1415
1416
      initElement(elInfo);
    }

1417
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1418
1419
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
1420
1421
1422
1423
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1424
      if (rememberElVec) {
1425
1426
1427
1428
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
1429
    
1430
1431
    ElementVector *vec = rememberElVec ? elementVector : userVec;

Thomas Witkowski's avatar
Thomas Witkowski committed
1432
    if (operat->uhOld && remember) {
1433
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
1434
      if (rememberElVec) {
1435
1436
1437
1438
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 
1439
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1440
    if (firstOrderAssemblerGrdPsi)
1441
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
1442
    if (zeroOrderAssembler)
1443
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
1444
        
Thomas Witkowski's avatar
Thomas Witkowski committed
1445
    if (rememberElVec) {
1446
      axpy(factor, *elementVector, *userVec);
1447
    }      
1448
1449
1450
1451
1452
1453
1454
1455
  }

  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
1456
1457
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
1458

1459
1460
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1461
    if (el != lastMatEl) {
1462
1463
      calculateElementMatrix(elInfo, NULL);
    }
1464
1465
1466
1467
1468
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
1469
1470
1471
      }
      (*vec)[i] += val;
    }
1472
1473
1474
    

    delete [] uhOldLoc;       
1475
1476
1477
1478
1479
1480
  }

  void Assembler::initElement(const ElInfo *elInfo, Quadrature *quad)
  {
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
1481
    if (secondOrderAssembler) 
1482
      secondOrderAssembler->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1483
    if (firstOrderAssemblerGrdPsi)
1484
      firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1485
    if (firstOrderAssemblerGrdPhi)
1486
      firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1487
    if (zeroOrderAssembler)
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
      zeroOrderAssembler->initElement(elInfo, quad);
  }

  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
					 const FiniteElemSpace *rowFESpace_,
					 const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    bool opt = (rowFESpace_ == colFESpace_);

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
1503
    secondOrderAssembler = 
1504
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
1505

Thomas Witkowski's avatar
Thomas Witkowski committed
1506
    firstOrderAssemblerGrdPsi = 
1507
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
1508

Thomas Witkowski's avatar
Thomas Witkowski committed
1509
    firstOrderAssemblerGrdPhi = 
1510
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
1511

Thomas Witkowski's avatar
Thomas Witkowski committed
1512
    zeroOrderAssembler = 
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
  }

  StandardAssembler::StandardAssembler(Operator  *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *rowFESpace_,
				       const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    remember = false;

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
1528
    secondOrderAssembler = 
1529
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1530
    firstOrderAssemblerGrdPsi = 
1531
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1532
    firstOrderAssemblerGrdPhi = 
1533
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1534
    zeroOrderAssembler = 
1535
1536
1537
1538
1539
1540
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

  ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1541
    if (!elMat) {
1542
1543
1544
1545
      elMat = NEW ElementMatrix(nRow, nCol);
    }

    elMat->set(0.0);
1546
      
1547
1548
    Element *element = elInfo->getElement();

1549
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1550
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element,
Thomas Witkowski's avatar
Thomas Witkowski committed
1551
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
1552
						   &(elMat->rowIndices));
1553

1554

Thomas Witkowski's avatar
Thomas Witkowski committed
1555
1556
1557
1558
1559
1560
1561
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
      colFESpace->getBasisFcts()->getLocalIndicesVec(element,
						     colFESpace->getAdmin(),
						     &(elMat->colIndices));
    }
1562
    
1563
1564
1565
1566
1567
1568
    return elMat;
  }

  ElementVector *Assembler::initElementVector(ElementVector *elVec, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1569
    if (!elVec) {
1570
1571
1572
1573
1574
1575
1576
1577
1578
      elVec = NEW ElementVector(nRow);
    }

    elVec->set(0.0);
  
    DOFAdmin *rowAdmin = rowFESpace->getAdmin();

    Element *element = elInfo->getElement();

1579
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
1580
1581
1582
1583
1584

    return elVec;
  }

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
1585
1586
  { 
    if (secondOrderAssembler) {
1587
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1588
1589
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1590
1591
1592
1593
1594
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1595
    if (firstOrderAssemblerGrdPsi) {
1596
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1597
1598
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1599
1600
1601
1602
1603
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1604
    if (firstOrderAssemblerGrdPhi) {
1605
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1606
1607
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1608
1609
1610
1611
1612
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1613
    if (zeroOrderAssembler) {
1614
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1615
1616
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1617
1618
1619
1620
1621
1622
1623
1624
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

}
For faster browsing, not all history is shown. View entire blame