Assembler.cc 46.1 KB
Newer Older
1001
1002
    int nPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,nPoints,NO_INIT);
1003
    int myRank = omp_get_thread_num();
1004

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

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

1017
      for (int i = 0; i < nRow; i++) {
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
	(*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;

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

    const int **nEntries = q01->getNumberEntries();
1045
    int myRank = omp_get_thread_num();
1046
    Lb[0].set(0.0);
1047
1048
1049

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

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

1054
1055
1056
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	l = q01->getLVec(i, j);
1057
	values = q01->getValVec(i, j);
1058
	double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1059
1060
1061
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][l[m]];
	}
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
	(*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
1074
    if (firstCall) {
1075
1076
1077
1078
1079
1080
1081
1082
1083
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int *nEntries = q1->getNumberEntries();
1084
    int myRank = omp_get_thread_num();
1085
    Lb[0].set(0.0);
1086
1087
1088

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

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

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

  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
  {
    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];
    }
  }
1124
1125
1126
1127
1128
1129
1130

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

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

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

1140
    tmpMat *= elInfo->getDet();
1141
1142
1143
    nEntries = q11->getNumberEntries();

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

	for (int j = i + 1; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
1157
	  values = q11->getValVec(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
1158
1159
	  val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
1160
	    val += values[m] * tmpMat[k[m]][l[m]];
Thomas Witkowski's avatar
Thomas Witkowski committed
1161
	  }
1162
1163
1164
1165
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1166
1167
1168
1169
1170
    } 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);
1171
	  values = q11->getValVec(i, j);
1172
	  double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1173
	  for (int m = 0; m < nEntries[i][j]; m++) {
1174
	    val += values[m] * tmpMat[k[m]][l[m]];
Thomas Witkowski's avatar
Thomas Witkowski committed
1175
	  }
1176
1177
1178
1179
1180
1181
1182
1183
	  (*mat)[i][j] += val;
	}
      }
    }
  }

  Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
1184
1185
1186
1187
1188
1189
  {
    tmpLALt.resize(omp_get_max_threads());
  }

  Quad2::~Quad2()
  {
1190
1191
1192
1193
1194
1195
1196
    if (!firstCall) {
      int nPoints = quadrature->getNumPoints();
      for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
	for (int j = 0; j < nPoints; j++) {
	  DELETE tmpLALt[i][j];
	}
	DELETE [] tmpLALt[i];
1197
1198
1199
      }
    }
  }
1200
1201
1202

  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
1203
1204
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
1205

Thomas Witkowski's avatar
Thomas Witkowski committed
1206
    if (firstCall) {
1207
1208
1209
1210
1211
      tmpLALt[myRank] = NEW DimMat<double>*[nPoints];
      for (int j = 0; j < nPoints; j++) {
	tmpLALt[myRank][j] = NEW DimMat<double>(dim, NO_INIT);
      }
    
1212
1213
1214
1215
1216
1217
1218
      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;
    }

1219
    DimMat<double> **LALt = tmpLALt[myRank];
1220
    
1221
1222
    for (int i = 0; i < nPoints; i++) {
      LALt[i]->set(0.0);
1223
    }
1224
    
1225
1226
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
1227
    }
1228
1229
        
    VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;
1230

1231
    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1232
      for (int iq = 0; iq < nPoints; iq++) {
1233
1234
1235
1236
1237
	(*LALt[iq]) *= elInfo->getDet();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1238
	for (int i = 0; i < nRow; i++) {
1239
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
1240
	    xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[i]);
1241

Thomas Witkowski's avatar
Thomas Witkowski committed
1242
	  for (int j = i + 1; j < nCol; j++) {
1243
1244
	    double val = quadrature->getWeight(iq) * 
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
1245
1246
1247
1248
1249
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
1250
    } else {      /*  non symmetric assembling   */
Thomas Witkowski's avatar
Thomas Witkowski committed
1251
      for (int iq = 0; iq < nPoints; iq++) {
1252
1253
1254
1255
1256
	(*LALt[iq]) *= elInfo->getDet();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1257
1258
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
1259
1260
	    (*mat)[i][j] += quadrature->getWeight(iq) * 
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
	  }
	}
      }
    }
  }

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

  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    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
1282
1283
    for (int iq = 0; iq < nPoints; iq++) {
      LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
1284
1285
      LALt[iq]->set(0.0);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1286

1287
1288
1289
1290
    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);
1291
1292
1293
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1297
	for (int i = 0; i < nCol; i++) {
1298
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
1299
1300
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
1301
	for (int i = 0; i < nRow; i++) {
1302
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
1303
1304
1305
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));

Thomas Witkowski's avatar
Thomas Witkowski committed
1306
	  for (int j = i + 1; j < nCol; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1307
	    double val = quadrature->getWeight(iq) * (grdPsi * ((*LALt[iq]) * grdPhi[j]));
1308
1309
1310
1311
1312
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
1313
1314
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
1315
1316
	(*LALt[iq]) *= elInfo->getDet();

Thomas Witkowski's avatar
Thomas Witkowski committed
1317
	for (int i = 0; i < nCol; i++) {
1318
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
1319
1320
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
1321
	for (int i = 0; i < nRow; i++) {
1322
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
Thomas Witkowski's avatar
Thomas Witkowski committed
1323
	  for (int j = 0; j < nCol; j++) {
1324
1325
1326
1327
1328
1329
1330
	    (*mat)[i][j] += quadrature->getWeight(iq) *
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	  }
	}
      }
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1331
1332
1333
    for (int iq = 0; iq < nPoints; iq++) 
      DELETE LALt[iq];

1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
    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
1354
  {}
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370

  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
1371
    
1372
1373
1374
    checkForNewTraverse();

    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
1375
    
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
    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
1403
    if (rememberElMat && userMat) {
1404
      axpy(factor, *elementMatrix, *userMat);
1405
    }      
1406
1407
1408
1409
1410
1411
1412
1413
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1414
    if (remember && factor != 1.0) {
1415
1416
1417
      rememberElVec = true;
    }

1418
    if (rememberElVec && !elementVector) {
1419
      elementVector = NEW ElementVector(nRow);
1420
    }
1421
1422
1423
1424
1425

    Element *el = elInfo->getElement();
    checkForNewTraverse();
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
1426
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
1427
1428
      initElement(elInfo);
    }
1429
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1430
1431
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
1432
1433
1434
1435
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1436
      if (rememberElVec) {
1437
1438
1439
1440
1441
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
    ElementVector *vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
1442
    if (operat->uhOld && remember) {
1443
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
1444
      if (rememberElVec) {
1445
1446
1447
1448
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 
1449
    if (firstOrderAssemblerGrdPsi) {
1450
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
1451
1452
    }
    if (zeroOrderAssembler) {
1453
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
1454
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1455
    if (rememberElVec) {
1456
      axpy(factor, *elementVector, *userVec);
1457
    }      
1458
1459
1460
1461
1462
1463
1464
1465
  }

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

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

1469
1470
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1471
    if (el != lastMatEl) {
1472
1473
      calculateElementMatrix(elInfo, NULL);
    }
1474
1475
1476
1477
1478
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
1479
1480
1481
      }
      (*vec)[i] += val;
    }
1482
1483
1484
    

    delete [] uhOldLoc;       
1485
1486
1487
1488
1489
1490
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1491
    if (secondOrderAssembler) 
1492
      secondOrderAssembler->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1493
    if (firstOrderAssemblerGrdPsi)
1494
      firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1495
    if (firstOrderAssemblerGrdPhi)
1496
      firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1497
    if (zeroOrderAssembler)
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
      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
1513
    secondOrderAssembler = 
1514
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
1515
    firstOrderAssemblerGrdPsi = 
1516
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
1517
    firstOrderAssemblerGrdPhi = 
1518
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
1519
    zeroOrderAssembler = 
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
      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
1535
    secondOrderAssembler = 
1536
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1537
    firstOrderAssemblerGrdPsi = 
1538
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1539
    firstOrderAssemblerGrdPhi = 
1540
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1541
    zeroOrderAssembler = 
1542
1543
1544
1545
1546
1547
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

  ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1548
    if (!elMat) {
1549
1550
1551
1552
      elMat = NEW ElementMatrix(nRow, nCol);
    }

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1553
1554
         
    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
1555
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
1556
						   &(elMat->rowIndices));
1557

1558

Thomas Witkowski's avatar
Thomas Witkowski committed
1559
1560
1561
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1562
      colFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
1563
1564
1565
						     colFESpace->getAdmin(),
						     &(elMat->colIndices));
    }
1566
    
1567
1568
1569
1570
1571
1572
    return elMat;
  }

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

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

    Element *element = elInfo->getElement();

1583
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
1584
1585
1586
1587
1588

    return elVec;
  }

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