Assembler.cc 44.2 KB
Newer Older
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
      }
    }
  }

  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;
    double val;

1018
    if (firstCall) {
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
      q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int **nEntries = q01->getNumberEntries();

    Lb[0].set(0.0);
1030
    for (int i = 0; i < static_cast<int>( terms.size()); i++) {
1031
1032
1033
1034
1035
      (static_cast<FirstOrderTerm*>((terms[i])))->getLb(elInfo, 1, Lb);
    }

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

1036
1037
1038
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	l = q01->getLVec(i, j);
1039
	values = q01->getValVec(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
1040
1041
1042
1043
	val = 0.0;
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][l[m]];
	}
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
	(*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;
    double val;

Thomas Witkowski's avatar
Thomas Witkowski committed
1057
    if (firstCall) {
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int *nEntries = q1->getNumberEntries();

    Lb[0].set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1069
    for (int i = 0; i < static_cast<int>(terms.size()); i++) {
1070
1071
1072
1073
1074
      (static_cast<FirstOrderTerm*>(terms[i]))->getLb(elInfo, 1, Lb);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1075
1076
    for (int i = 0; i < nRow; i++) {
      k = q1->getKVec(i);
1077
      values = q1->getValVec(i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1078
1079
1080
1081
      val = 0.0;
      for (int m = 0; m < nEntries[i]; m++) {
	val += values[m] * Lb[0][k[m]];
      }
1082
1083
1084
1085
1086
1087
      (*vec)[i] += val;
    }
  }

  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
Thomas Witkowski's avatar
Thomas Witkowski committed
1088
  {}
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098

  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    DimMat<double> **LALt = NEW DimMat<double>*;
    *LALt=NEW DimMat<double>(dim, NO_INIT);
    const int **nEntries;
    const int *k, *l;
    const double *values;
    double val;

Thomas Witkowski's avatar
Thomas Witkowski committed
1099
    if (firstCall) {
1100
1101
1102
1103
1104
1105
1106
1107
      q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      firstCall = false;
    }

    LALt[0]->set(0.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
1108
    for (int i = 0; i < static_cast<int>( terms.size()); i++) {
1109
1110
1111
1112
1113
1114
1115
1116
      (static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, 1, LALt);
    }

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

    nEntries = q11->getNumberEntries();

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1117
1118
1119
      for (int i = 0; i < nRow; i++) {
	k = q11->getKVec(i, i);
	l = q11->getLVec(i, i);
1120
	values = q11->getValVec(i, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1121
1122
1123
1124
	val = 0.0;
	for (int m = 0; m < nEntries[i][i]; m++) {
	  val += values[m] * (*LALt[0])[k[m]][l[m]];
	}
1125
	(*mat)[i][i] += val;
Thomas Witkowski's avatar
Thomas Witkowski committed
1126
1127
1128
1129

	for (int j = i + 1; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
1130
	  values = q11->getValVec(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
1131
1132
1133
1134
	  val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
	    val += values[m] * (*LALt[0])[k[m]][l[m]];
	  }
1135
1136
1137
1138
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1139
1140
1141
1142
1143
    } 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);
1144
	  values = q11->getValVec(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
1145
1146
1147
1148
	  val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
	    val += values[m] * (*LALt[0])[k[m]][l[m]];
	  }
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
	  (*mat)[i][j] += val;
	}
      }
    }

    DELETE *LALt;
    DELETE LALt;
  }

  Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
Thomas Witkowski's avatar
Thomas Witkowski committed
1160
  {}
1161
1162
1163
1164
1165
1166

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1167
    if (firstCall) {
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
      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];
Thomas Witkowski's avatar
Thomas Witkowski committed
1178
1179
1180
1181
    for (int i = 0; i < nPoints;i++) {
      LALt[i] = NEW DimMat<double>(dim, NO_INIT);
    }
    for (int iq = 0; iq < nPoints; iq++) {
1182
1183
      LALt[iq]->set(0.0);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1184
    for (int i = 0; i < static_cast<int>(terms.size()); i++) {
1185
1186
1187
1188
      (static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, nPoints, LALt);
    }

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1189
      for (int iq = 0; iq < nPoints; iq++) {
1190
1191
1192
1193
1194
	(*LALt[iq]) *= elInfo->getDet();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1195
	for (int i = 0; i < nRow; i++) {
1196
1197
1198
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
	    ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));

Thomas Witkowski's avatar
Thomas Witkowski committed
1199
	  for (int j = i + 1; j < nCol; j++) {
1200
1201
1202
1203
1204
1205
1206
1207
	    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
1208
      for (int iq = 0; iq < nPoints; iq++) {
1209
1210
1211
1212
1213
	(*LALt[iq]) *= elInfo->getDet();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1214
1215
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
1216
1217
1218
1219
1220
1221
1222
	    (*mat)[i][j] += quadrature->getWeight(iq) *
	      ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[j]));
	  }
	}
      }
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1223
1224
1225
    for (int i = 0;i < nPoints; i++) 
      DELETE LALt[i];

1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
    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
1245
1246
    for (int iq = 0; iq < nPoints; iq++) {
      LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
1247
1248
      LALt[iq]->set(0.0);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1249
1250

    for (int i = 0; i < static_cast<int>(terms.size()); i++) {
1251
1252
1253
1254
      (static_cast<SecondOrderTerm*>(terms[i]))->getLALt(elInfo, nPoints, LALt);
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1258
	for (int i = 0; i < nCol; i++) {
1259
1260
1261
	  grdPhi[i] = (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq));
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
1262
	for (int i = 0; i < nRow; i++) {
1263
1264
1265
1266
	  grdPsi = (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq));
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));

Thomas Witkowski's avatar
Thomas Witkowski committed
1267
	  for (int j = i + 1; j < nCol; j++) {
1268
1269
1270
1271
1272
1273
	    val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1274
1275
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
1276
1277
	(*LALt[iq]) *= elInfo->getDet();

Thomas Witkowski's avatar
Thomas Witkowski committed
1278
	for (int i = 0; i < nCol; i++) {
1279
1280
1281
	  grdPhi[i] = (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq));
	}

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1292
1293
1294
    for (int iq = 0; iq < nPoints; iq++) 
      DELETE LALt[iq];

1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
    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
1315
  {}
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331

  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
1332
    
1333
1334
1335
    checkForNewTraverse();

    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
1336
    
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
    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
1364
    if (rememberElMat && userMat) {
1365
      axpy(factor, *elementMatrix, *userMat);
Thomas Witkowski's avatar
Thomas Witkowski committed
1366
    }    
1367
1368
1369
1370
1371
1372
1373
1374
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1375
    if (remember && factor != 1.0) {
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
      rememberElVec = true;
    }

    if(rememberElVec && !elementVector)
      elementVector = NEW ElementVector(nRow);

    Element *el = elInfo->getElement();

    checkForNewTraverse();

    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
1388
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
1389
1390
1391
      initElement(elInfo);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1392
1393
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
1394
1395
1396
1397
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1398
      if (rememberElVec) {
1399
1400
1401
1402
1403
1404
1405
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }

    ElementVector *vec = rememberElVec ? elementVector : userVec;

Thomas Witkowski's avatar
Thomas Witkowski committed
1406
    if (operat->uhOld && remember) {
1407
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
1408
      if (rememberElVec) {
1409
1410
1411
1412
1413
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 

Thomas Witkowski's avatar
Thomas Witkowski committed
1414
    if (firstOrderAssemblerGrdPsi)
1415
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
1416
    if (zeroOrderAssembler)
1417
1418
      zeroOrderAssembler->calculateElementVector(elInfo, vec);

Thomas Witkowski's avatar
Thomas Witkowski committed
1419
    if (rememberElVec) {
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
      axpy(factor, *elementVector, *userVec);
    }  
  }

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

    Element *el = elInfo->getElement(); 
    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
    const double *uhOldLoc = operat->uhOld->getLocalVector(el, NULL);
    int n = basFcts->getNumber();

Thomas Witkowski's avatar
Thomas Witkowski committed
1433
    if (el != lastMatEl) {
1434
1435
1436
1437
      calculateElementMatrix(elInfo, NULL);
    }

    double val;
Thomas Witkowski's avatar
Thomas Witkowski committed
1438
1439
1440
    for (int i = 0; i < n; i++) {
      val = 0.0;
      for (int j = 0; j < n; j++) {
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
	val += (*elementMatrix)[i][j]*uhOldLoc[j];
      }
      (*vec)[i] += val;
    }
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1451
    if (secondOrderAssembler) 
1452
      secondOrderAssembler->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1453
    if (firstOrderAssemblerGrdPsi)
1454
      firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1455
    if (firstOrderAssemblerGrdPhi)
1456
      firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1457
    if (zeroOrderAssembler)
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
      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
1473
    secondOrderAssembler = 
1474
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
1475
    firstOrderAssemblerGrdPsi = 
1476
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
1477
    firstOrderAssemblerGrdPhi = 
1478
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
1479
    zeroOrderAssembler = 
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
      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
1495
    secondOrderAssembler = 
1496
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1497
    firstOrderAssemblerGrdPsi = 
1498
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1499
    firstOrderAssemblerGrdPhi = 
1500
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1501
    zeroOrderAssembler = 
1502
1503
1504
1505
1506
1507
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

  ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1508
    if (!elMat) {
1509
1510
1511
1512
      elMat = NEW ElementMatrix(nRow, nCol);
    }

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1513
        
1514
1515
    Element *element = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
1516
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element,
Thomas Witkowski's avatar
Thomas Witkowski committed
1517
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
1518
						   &(elMat->rowIndices));
1519

Thomas Witkowski's avatar
Thomas Witkowski committed
1520
1521
1522
1523
1524
1525
1526
1527
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
      colFESpace->getBasisFcts()->getLocalIndicesVec(element,
						     colFESpace->getAdmin(),
						     &(elMat->colIndices));
    }

1528
1529
1530
1531
1532
1533
    return elMat;
  }

  ElementVector *Assembler::initElementVector(ElementVector *elVec, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1534
    if (!elVec) {
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
      elVec = NEW ElementVector(nRow);
    }

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

    Element *element = elInfo->getElement();

    elVec->dofIndices = 
      rowFESpace->getBasisFcts()->getLocalIndices(element,
						  rowAdmin,
						  NULL);

    return elVec;
  }

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
1553
1554
  { 
    if (secondOrderAssembler) {
1555
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1556
1557
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1558
1559
1560
1561
1562
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1563
    if (firstOrderAssemblerGrdPsi) {
1564
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1565
1566
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1567
1568
1569
1570
1571
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1572
    if (firstOrderAssemblerGrdPhi) {
1573
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1574
1575
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1576
1577
1578
1579
1580
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1581
    if (zeroOrderAssembler) {
1582
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
1583
1584
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
1585
1586
1587
1588
1589
1590
1591
1592
	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