Assembler.cc 45.5 KB
Newer Older
1001
    int myRank = omp_get_thread_num();
1002

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

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

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

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

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

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

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

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

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

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

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

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

  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
Thomas Witkowski's avatar
Thomas Witkowski committed
1104
  {}
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114

  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
1115
    if (firstCall) {
1116
1117
1118
1119
1120
1121
1122
      q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      firstCall = false;
    }

    LALt[0]->set(0.0);
1123
    int myRank = omp_get_thread_num();
1124

1125
1126
    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
1127
1128
1129
1130
1131
1132
1133
    }

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

    nEntries = q11->getNumberEntries();

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1134
1135
1136
      for (int i = 0; i < nRow; i++) {
	k = q11->getKVec(i, i);
	l = q11->getLVec(i, i);
1137
	values = q11->getValVec(i, i);
Thomas Witkowski's avatar
Thomas Witkowski committed
1138
1139
1140
1141
	val = 0.0;
	for (int m = 0; m < nEntries[i][i]; m++) {
	  val += values[m] * (*LALt[0])[k[m]][l[m]];
	}
1142
	(*mat)[i][i] += val;
Thomas Witkowski's avatar
Thomas Witkowski committed
1143
1144
1145
1146

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1184
    if (firstCall) {
1185
1186
1187
1188
1189
1190
1191
1192
1193
      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];
1194
1195
    int myRank = omp_get_thread_num();

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

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1207
      for (int iq = 0; iq < nPoints; iq++) {
1208
1209
1210
1211
1212
	(*LALt[iq]) *= elInfo->getDet();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1213
	for (int i = 0; i < nRow; i++) {
1214
1215
1216
	  (*mat)[i][i] += quadrature->getWeight(iq) * 
	    ((*grdPsi)[i] * ((*LALt[iq]) * (*grdPhi)[i]));

Thomas Witkowski's avatar
Thomas Witkowski committed
1217
	  for (int j = i + 1; j < nCol; j++) {
1218
1219
1220
1221
1222
1223
1224
1225
	    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
1226
      for (int iq = 0; iq < nPoints; iq++) {
1227
1228
1229
1230
1231
	(*LALt[iq]) *= elInfo->getDet();

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

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

1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
    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
1263
1264
    for (int iq = 0; iq < nPoints; iq++) {
      LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
1265
1266
      LALt[iq]->set(0.0);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1267

1268
1269
1270
1271
    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);
1272
1273
1274
    }

    if (symmetric) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1275
      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
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
1280
1281
	}

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1312
1313
1314
    for (int iq = 0; iq < nPoints; iq++) 
      DELETE LALt[iq];

1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
    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
1335
  {}
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351

  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
1352
    
1353
1354
1355
    checkForNewTraverse();

    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
1356
    
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
    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
1384
    if (rememberElMat && userMat) {
1385
      axpy(factor, *elementMatrix, *userMat);
1386
    }      
1387
1388
1389
1390
1391
1392
1393
1394
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1395
    if (remember && factor != 1.0) {
1396
1397
1398
      rememberElVec = true;
    }

1399
    if (rememberElVec && !elementVector)
1400
1401
1402
1403
1404
1405
1406
1407
      elementVector = NEW ElementVector(nRow);

    Element *el = elInfo->getElement();

    checkForNewTraverse();

    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
1408
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
1409
1410
1411
      initElement(elInfo);
    }

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

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

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

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

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

    delete [] uhOldLoc;       
1470
1471
1472
1473
1474
1475
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1476
    if (secondOrderAssembler) 
1477
      secondOrderAssembler->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1478
    if (firstOrderAssemblerGrdPsi)
1479
      firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1480
    if (firstOrderAssemblerGrdPhi)
1481
      firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
1482
    if (zeroOrderAssembler)
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
      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
1498
    secondOrderAssembler = 
1499
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
1500

Thomas Witkowski's avatar
Thomas Witkowski committed
1501
    firstOrderAssemblerGrdPsi = 
1502
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
1503

Thomas Witkowski's avatar
Thomas Witkowski committed
1504
    firstOrderAssemblerGrdPhi = 
1505
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
1506

Thomas Witkowski's avatar
Thomas Witkowski committed
1507
    zeroOrderAssembler = 
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
      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
1523
    secondOrderAssembler = 
1524
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1525
    firstOrderAssemblerGrdPsi = 
1526
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1527
    firstOrderAssemblerGrdPhi = 
1528
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
1529
    zeroOrderAssembler = 
1530
1531
1532
1533
1534
1535
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

  ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1536
    if (!elMat) {
1537
1538
1539
1540
      elMat = NEW ElementMatrix(nRow, nCol);
    }

    elMat->set(0.0);
1541
      
1542
1543
    Element *element = elInfo->getElement();

1544
    
Thomas Witkowski's avatar
Thomas Witkowski committed
1545
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element,
Thomas Witkowski's avatar
Thomas Witkowski committed
1546
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
1547
						   &(elMat->rowIndices));
1548

1549

Thomas Witkowski's avatar
Thomas Witkowski committed
1550
1551
1552
1553
1554
1555
1556
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
      colFESpace->getBasisFcts()->getLocalIndicesVec(element,
						     colFESpace->getAdmin(),
						     &(elMat->colIndices));
    }
1557
    
1558
1559
1560
1561
1562
1563
    return elMat;
  }

  ElementVector *Assembler::initElementVector(ElementVector *elVec, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1564
    if (!elVec) {
1565
1566
1567
1568
1569
1570
1571
1572
1573
      elVec = NEW ElementVector(nRow);
    }

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

    Element *element = elInfo->getElement();

1574
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
1575
1576
1577
1578
1579

    return elVec;
  }

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