DOFVector.hh 57.3 KB
Newer Older
1001

1002

1003
1004
1005
1006
  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
1007
1008
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
Thomas Witkowski's avatar
Thomas Witkowski committed
1009
	 op != this->operators.end(); ++op)
1010
      fillFlag |= (*op)->getFillFlag();
Thomas Witkowski's avatar
Thomas Witkowski committed
1011

1012
1013
1014
    return fillFlag;
  }

1015

Thomas Witkowski's avatar
Thomas Witkowski committed
1016
1017
1018
1019
1020
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
1021
	 it != this->operators.end(); ++it)
Thomas Witkowski's avatar
Thomas Witkowski committed
1022
1023
1024
      (*it)->finishAssembling();
  }

1025

1026
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
1027
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs)
1028
  {
1029
    this->feSpace = rhs.feSpace;
1030
    this->dim = rhs.dim;
1031
    this->nBasFcts = rhs.nBasFcts;
1032
    vec = rhs.vec;
1033
    this->elementVector.change_dim(this->nBasFcts);
1034
    coarsenOperation = rhs.coarsenOperation;
1035
    refineOperation = rhs.refineOperation;
1036
1037
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
1038
    
1039
    if (rhs.boundaryManager) {
1040
1041
1042
      if (this->boundaryManager) 
	delete this->boundaryManager; 

1043
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
1044
    } else {
1045
      this->boundaryManager = nullptr;
1046
1047
    }

1048
1049
1050
    return *this;
  }

1051

1052
1053
1054
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
1055
    FUNCNAME_DBG("DOFVector<T>::operator*=(DOFVector<T>& x, T scal)");
1056

1057
    TEST_EXIT_DBG(x.getFeSpace() && x.getFeSpace()->getAdmin())
1058
      ("pointer is nullptr: %8X, %8X\n", x.getFeSpace(), x.getFeSpace()->getAdmin());
1059

Thomas Witkowski's avatar
Thomas Witkowski committed
1060
1061
    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), 
						USED_DOFS);
1062
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
1063
      (*vecIterator) *= scal; 
1064

1065
1066
1067
1068
1069
1070
1071
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
1072
    FUNCNAME_DBG("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
1073
    
1074
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
1075
      ("feSpace is nullptr: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
1076
1077
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
1078
      ("no admin or different admins: %8X, %8X\n",
1079
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1080
1081
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
1082
1083
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1084
1085
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
1086
1087
      *xIterator += *yIterator; 

1088
1089
1090
    return x;
  }

1091

1092
1093
1094
  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y)
  {
1095
    FUNCNAME_DBG("DOFVector<T>::operator-=(DOFVector<T>& x, const DOFVector<T>& y)");
1096

1097
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
1098
      ("feSpace is nullptr: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
1099
1100
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
1101
      ("no admin or different admins: %8X, %8X\n",
1102
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1103
1104
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
1105
1106
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1107
1108
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
1109
      *xIterator -= *yIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
1110

1111
1112
1113
    return x;
  }

1114

1115
1116
1117
  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
1118
    FUNCNAME_DBG("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
1119
    
1120
    TEST_EXIT_DBG(x.getFeSpace() && y.getFeSpace())
1121
      ("feSpace is nullptr: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
1122
1123
    TEST_EXIT_DBG(x.getFeSpace()->getAdmin() &&
	      (x.getFeSpace()->getAdmin() == y.getFeSpace()->getAdmin()))
1124
      ("no admin or different admins: %8X, %8X\n",
1125
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1126
1127
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
1128
1129
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1130
1131
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
	 ++xIterator, ++yIterator)
1132
      *xIterator *= *yIterator; 
1133
1134
1135
1136

    return x;
  }

1137

1138
1139
1140
1141
  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");
1142
    const DOFAdmin *admin = nullptr;
1143

1144
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
1145
      ("feSpace is nullptr: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
1146
    TEST_EXIT((admin = x.getFeSpace()->getAdmin()) && (admin == y.getFeSpace()->getAdmin()))
1147
      ("no admin or different admins: %8X, %8X\n",
1148
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1149
1150
    TEST_EXIT(x.getSize() == y.getSize())("different sizes\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
1151
    T dot = 0;
1152
1153
1154

    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1155
    for (xIterator.reset(), yIterator.reset(); !xIterator.end();
1156
	 ++xIterator, ++yIterator) 
Thomas Witkowski's avatar
Thomas Witkowski committed
1157
      dot += (*xIterator) * (*yIterator);
1158

Thomas Witkowski's avatar
Thomas Witkowski committed
1159
    return dot;
1160
1161
  }

1162

1163
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
1164
1165
  void mv(MatrixTranspose transpose, const DOFMatrix &a, const DOFVector<T>&x,
	  DOFVector<T> &result, bool add)
1166
  {
1167
1168
1169
    // Unfortunately needed
    using namespace mtl;

1170
1171
    FUNCNAME("DOFVector<T>::mv()");

1172
1173
    TEST_EXIT(a.getRowFeSpace() && a.getColFeSpace() && 
	      x.getFeSpace() && result.getFeSpace())
1174
      ("getFeSpace() is nullptr: %8X, %8X, %8X, %8X\n", 
1175
       a.getRowFeSpace(), a.getColFeSpace(), x.getFeSpace(), result.getFeSpace());
1176

1177
1178
    const DOFAdmin *rowAdmin = a.getRowFeSpace()->getAdmin();
    const DOFAdmin *colAdmin = a.getColFeSpace()->getAdmin();
1179
1180

    TEST_EXIT((rowAdmin && colAdmin) && 
1181
1182
1183
1184
	      (((transpose == NoTranspose) && (colAdmin == x.getFeSpace()->getAdmin()) && 
		(rowAdmin == result.getFeSpace()->getAdmin()))||
	       ((transpose == Transpose) && (rowAdmin == x.getFeSpace()->getAdmin()) && 
		(colAdmin == result.getFeSpace()->getAdmin()))))
1185
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
1186
1187
       a.getRowFeSpace()->getAdmin(), a.getColFeSpace()->getAdmin(), 
       x.getFeSpace()->getAdmin(), result.getFeSpace()->getAdmin());
1188
1189


1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
    if (transpose == NoTranspose) {
      mtl::dense_vector<T> tmp_x(x.getUsedSize());
      mtl::dense_vector<T> tmp_result(result.getUsedSize());

      {
	int counter = 0;
	typename DOFVector<T>::Iterator it(const_cast<DOFVector<T>*>(&x), USED_DOFS);
	for (it.reset(); !it.end(); ++it)
	  tmp_x[counter++] = *it;
      }

      mult(a.getBaseMatrix(), tmp_x, tmp_result);

      {
	int counter = 0;
	typename DOFVector<T>::Iterator it(&result, USED_DOFS);
	for (it.reset(); !it.end(); ++it)
	  if (add)
	    *it += tmp_result[counter++];
	  else
	    *it = tmp_result[counter++];
      }


    } else if (transpose == Transpose) {
      ERROR_EXIT("Not yet implemented!\n");

//       mult(trans(const_cast<DOFMatrix::base_matrix_type&>(a.getBaseMatrix())),
// 	   x, result);
    }
1220
    else 
1221
1222
1223
      ERROR_EXIT("transpose = %d\n", transpose);
  }

1224

1225
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
1226
  void axpy(double alpha, const DOFVector<T>& x, DOFVector<T>& y)
1227
1228
1229
  {
    FUNCNAME("DOFVector<T>::axpy()");

1230
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
1231
      ("feSpace is nullptr: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
1232

1233
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
1234

1235
    TEST_EXIT((admin) && (admin == y.getFeSpace()->getAdmin()))
1236
      ("no admin or different admins: %8X, %8X\n",
1237
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1238
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
1239
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
1240
1241
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
1242
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
1243
1244
       admin->getUsedSize());

1245
1246
1247
    for (int i = 0; i < y.getSize(); i++)
      if (!admin->isDofFree(i))
	y[i] += alpha * x[i];
1248
1249
  }

1250

1251
1252
1253
  template<typename T>
  const DOFVector<T>& operator*(const DOFVector<T>& v, double d)
  {
1254
1255
    static DOFVector<T> result; // TODO: REMOVE STATIC
    return mult(d, v, result); 
1256
1257
  }

1258

1259
1260
1261
  template<typename T>
  const DOFVector<T>& operator*(double d, const DOFVector<T>& v)
  {
1262
    static DOFVector<T> result; // TODO: REMOVE STATIC
1263
1264
1265
    return mult(d, v, result);
  }

1266

1267
1268
1269
  template<typename T>
  const DOFVector<T>& operator+(const DOFVector<T> &v1 , const DOFVector<T> &v2)
  {
1270
    static DOFVector<T> result; // TODO: REMOVE STATIC
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
    return add(v1, v2, result);
  }


  template<typename T>
  void xpay(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::xpay()");

1282
    TEST_EXIT(x.getFeSpace() && y.getFeSpace())
1283
      ("feSpace is nullptr: %8X, %8X\n", x.getFeSpace(), y.getFeSpace());
1284

1285
    const DOFAdmin *admin = x.getFeSpace()->getAdmin();
1286

1287
    TEST_EXIT(admin && (admin == y.getFeSpace()->getAdmin()))
1288
      ("no admin or different admins: %8X, %8X\n",
1289
       x.getFeSpace()->getAdmin(), y.getFeSpace()->getAdmin());
1290
    TEST_EXIT(static_cast<int>(x.getSize()) >= admin->getUsedSize())
1291
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
1292
1293
       admin->getUsedSize());
    TEST_EXIT(static_cast<int>(y.getSize()) >= admin->getUsedSize())
1294
      ("y.size = %d too small: admin->size = %d\n", y.getSize(), 
1295
1296
       admin->getUsedSize());

1297
1298
1299
    for (int i = 0; i < y.getSize(); i++)
      if (!admin->isDofFree(i))
	y[i] = alpha * y[i] + x[i];
1300
1301
  }

1302

1303
1304
1305
1306
1307
1308
1309
  template<typename T>
  inline const DOFVector<T>& mult(double scal,
				  const DOFVector<T>& v,
				  DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1310
1311
1312
    for (vIterator.reset(), rIterator.reset();
	 !vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = scal * (*vIterator);
1313
1314
1315
1316

    return result;
  }

1317

1318
1319
1320
1321
1322
1323
1324
  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v,
				 double scal,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1325
1326
1327
1328
    for (vIterator.reset(), rIterator.reset();
	!vIterator.end(); ++vIterator, ++rIterator) 
      *rIterator = (*vIterator) + scal;

1329
1330
1331
    return result;
  }

1332

1333
  template<typename T>
1334
  inline const DOFVector<T>& add(const DOFVector<T>& v1, 
1335
1336
1337
1338
1339
1340
				 const DOFVector<T>& v2,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
    typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
1341
1342
1343
    for (v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	 !v1Iterator.end(); ++v1Iterator, ++v2Iterator, ++rIterator) 
      *rIterator = (*v1Iterator) + (*v2Iterator);
1344

Thomas Witkowski's avatar
Thomas Witkowski committed
1345
    return result;
1346
1347
  }

1348

1349
  template<typename T>
1350
1351
  void DOFVectorBase<T>::getLocalVector(const Element *el, 
					mtl::dense_vector<T>& d) const
1352
  {
1353
    FUNCNAME_DBG("DOFVectorBase<T>::getLocalVector()");
1354
1355
1356
1357

    TEST_EXIT_DBG(feSpace->getMesh() == el->getMesh())
      ("Element is defined on a different mesh than the DOF vector!\n");

1358
    std::vector<DegreeOfFreedom> localIndices(nBasFcts);
1359
    const DOFAdmin* admin = feSpace->getAdmin();
Thomas Witkowski's avatar
Blub    
Thomas Witkowski committed
1360
    feSpace->getBasisFcts()->getLocalIndices(el, admin, localIndices);
1361

1362
    for (int i = 0; i < nBasFcts; i++)
1363
      d[i] = (*this)[localIndices[i]];
1364
1365
  }

1366

1367
  template<typename T>
1368
1369
1370
1371
  void DOFVectorBase<T>::getVecAtQPs(const ElInfo *elInfo, 
				     const Quadrature *quad,
				     const FastQuadrature *quadFast,
				     mtl::dense_vector<T>& vecAtQPs) const
1372
  {
1373
    FUNCNAME_DBG("DOFVector<T>::getVecAtQPs()");
1374
1375
1376
1377
1378
    
    TEST_EXIT_DBG(quad || quadFast)
      ("Neither quad nor quadFast defined!\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");    
1379
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1380
      ("Invalid basis functions!");
1381

1382
1383
1384
1385
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    mtl::dense_vector<T> localVec(nBasFcts);
    getLocalVector(elInfo->getElement(), localVec);
1386
1387
    
    if (quadFast) {
1388
      // using precalculated values at QPs
1389
      const mtl::dense2D<double>& phi = quadFast->getPhi();
1390
1391
      vecAtQPs.change_dim(num_rows(phi)); 	// = quadrature->getNumPoints()
      vecAtQPs = phi * localVec; 		// Matrix<double> * Vector<T>
1392
    } else {
1393
1394
      // evaluate basisFunctions at QPs
      int nPoints = quad->getNumPoints();
1395
1396
      vecAtQPs.change_dim(nPoints);

1397
1398
      for (int iq = 0; iq < nPoints; iq++) {
       	vecAtQPs[iq] = 0.0;
1399
       	for (int j = 0; j < nBasFcts; j++)
1400
1401
       	  vecAtQPs[iq] += 
       	    localVec[j] * (*(basFcts->getPhi(j)))(quad->getLambda(iq));
1402
      }
1403
1404
1405
    }
  }

1406

1407
1408
1409
1410
1411
1412
1413
  template<typename T>
  void DOFVectorBase<T>::getVecAtQPs(const ElInfo *smallElInfo, 
					  const ElInfo *largeElInfo,
					  const Quadrature *quad,
					  const FastQuadrature *quadFast,
					  mtl::dense_vector<T>& vecAtQPs) const
  {
1414
    FUNCNAME_DBG("DOFVector<T>::getVecAtQPs()");
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
 
    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    if (smallElInfo->getMesh() == feSpace->getMesh())
      return getVecAtQPs(smallElInfo, quad, quadFast, vecAtQPs);    

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int nPoints = quad->getNumPoints();

    vecAtQPs.change_dim(nPoints);
    
    mtl::dense_vector<T> localVec(nBasFcts);
1432
    this->getLocalVector(largeElInfo->getElement(), localVec);
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());

    for (int iq = 0; iq < nPoints; iq++) {
      vecAtQPs[iq] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int k = 0; k < nBasFcts; k++)
	  val += m[j][k] * (*(basFcts->getPhi(k)))(quad->getLambda(iq));

	vecAtQPs[iq] += localVec[j] * val;
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getGrdAtQPs(const ElInfo *elInfo,
				     const Quadrature *quad,
				     const FastQuadrature *quadFast,
				     mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const
  {
1454
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int dow = Global::getGeo(WORLD);
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
1468
    this->getLocalVector(elInfo->getElement(), localVec);
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478

    mtl::dense_vector<T> grd1(dim + 1);
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

    grdAtQPs.change_dim(nPoints);
    if (quadFast) {
      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

Thomas Witkowski's avatar
Thomas Witkowski committed
1479
	for (int j = 0; j < nBasFcts; j++) { // #BasisFunctions 
1480
1481
	  for (int k = 0; k < parts; k++)  // #edges (2d) or #faces (3d)
	    grd1[k] += quadFast->getGradient(iq, j, k) * localVec[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
1482
	}
1483

Thomas Witkowski's avatar
Thomas Witkowski committed
1484
	for (int l = 0; l < dow; l++) {
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
	  nullify(grdAtQPs[iq][l]);
	  for (int k = 0; k < parts; k++)
	    grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
	}
      }

    } else {
      mtl::dense_vector<double> grdPhi(dim + 1);

      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) {
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(iq), grdPhi);
	  for (int k = 0; k < parts; k++)
	    grd1[k] += grdPhi[k] * localVec[j];
	}

	for (int l = 0; l < dow; l++) {
	  nullify(grdAtQPs[iq][l]);
	  for (int k = 0; k < parts; k++)
	    grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
	}
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getGrdAtQPs(const ElInfo *smallElInfo,
				     const ElInfo *largeElInfo,
				     const Quadrature *quad,
				     const FastQuadrature *quadFast,
				     mtl::dense_vector<typename GradientType<T>::type> &grdAtQPs) const
  {
1520
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    if (smallElInfo->getMesh() == feSpace->getMesh())
      return getGrdAtQPs(smallElInfo, quad, quadFast, grdAtQPs);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int dow = Global::getGeo(WORLD);
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
1537
    this->getLocalVector(largeElInfo->getElement(), localVec);
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580

    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());

    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    mtl::dense_vector<T> grd1(dim + 1);
    mtl::dense_vector<double> grdPhi(dim + 1);
    mtl::dense_vector<double> tmp(dim + 1);

    grdAtQPs.change_dim(nPoints);
    for (int iq = 0; iq < nPoints; iq++) {
      nullify(grd1);

      for (int j = 0; j < nBasFcts; j++) {
	grdPhi = 0.0;

	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(iq), tmp);
	  tmp *= m[j][k];
	  grdPhi += tmp;
	}

	for (int k = 0; k < parts; k++)
	  grd1[k] += grdPhi[k] * localVec[j];
      }

      for (int l = 0; l < dow; l++) {
	nullify(grdAtQPs[iq][l]);
	for (int k = 0; k < parts; k++)
	  grdAtQPs[iq][l] += grdLambda[k][l] * grd1[k];
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getDerivativeAtQPs(const ElInfo *elInfo,
					    const Quadrature *quad,
					    const FastQuadrature *quadFast,
					    int comp,
					    mtl::dense_vector<T> &derivativeAtQPs) const
  {
1581
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
1594
    this->getLocalVector(elInfo->getElement(), localVec);
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641

    mtl::dense_vector<T> grd1(dim + 1);
    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();

    derivativeAtQPs.change_dim(nPoints);
    if (quadFast) {
      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) // #BasisFunctions
	  for (int k = 0; k < parts; k++)  // #edges (2d) or #faces (3d)
	    grd1[k] += quadFast->getGradient(iq, j, k) * localVec[j];

	  nullify(derivativeAtQPs[iq]);
	  for (int k = 0; k < parts; k++)
	    derivativeAtQPs[iq] += grdLambda[k][comp] * grd1[k];
      }

    } else {
      mtl::dense_vector<double> grdPhi(dim + 1);

      for (int iq = 0; iq < nPoints; iq++) {
	nullify(grd1);

	for (int j = 0; j < nBasFcts; j++) {
	  (*(basFcts->getGrdPhi(j)))(quad->getLambda(iq), grdPhi);
	  for (int k = 0; k < parts; k++)
	    grd1[k] += grdPhi[k] * localVec[j];
	}

	  nullify(derivativeAtQPs[iq]);
	  for (int k = 0; k < parts; k++)
	    derivativeAtQPs[iq] += grdLambda[k][comp] * grd1[k];
      }
    }
  }


  template<typename T>
  void DOFVectorBase<T>::getDerivativeAtQPs(const ElInfo *smallElInfo,
					    const ElInfo *largeElInfo,
					    const Quadrature *quad,
					    const FastQuadrature *quadFast,
					    int comp,
					    mtl::dense_vector<T> &derivativeAtQPs) const
  {
1642
    FUNCNAME_DBG("DOFVector<T>::getGrdAtQPs()");
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657

    TEST_EXIT_DBG(quad || quadFast)("neither quad nor quadFast defined\n");
    TEST_EXIT_DBG(!(quad && quadFast) || quad == quadFast->getQuadrature())
      ("quad != quadFast->quadrature\n");
    TEST_EXIT_DBG(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
      ("invalid basis functions");

    if (smallElInfo->getMesh() == feSpace->getMesh())
      return getDerivativeAtQPs(smallElInfo, quad, quadFast, comp, derivativeAtQPs);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasFcts  = basFcts->getNumber();
    int nPoints = quadFast ? quadFast->getQuadrature()->getNumPoints() : quad->getNumPoints();

    mtl::dense_vector<T> localVec(nBasFcts);
1658
    this->getLocalVector(largeElInfo->getElement(), localVec);
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690

    mtl::dense2D<double> &m = smallElInfo->getSubElemCoordsMat(basFcts->getDegree());

    int parts = Global::getGeo(PARTS, dim);
    const DimVec<WorldVector<double> > &grdLambda = largeElInfo->getGrdLambda();

    mtl::dense_vector<T> grd1(dim + 1);
    mtl::dense_vector<double> grdPhi(dim + 1);
    mtl::dense_vector<double> tmp(dim + 1);

    derivativeAtQPs.change_dim(nPoints);
    for (int iq = 0; iq < nPoints; iq++) {
      nullify(grd1);

      for (int j = 0; j < nBasFcts; j++) {
	grdPhi = 0.0;

	for (int k = 0; k < nBasFcts; k++) {
	  (*(basFcts->getGrdPhi(k)))(quad->getLambda(iq), tmp);
	  tmp *= m[j][k];
	  grdPhi += tmp;
	}

	for (int k = 0; k < parts; k++)
	  grd1[k] += grdPhi[k] * localVec[j];
      }

      nullify(derivativeAtQPs[iq]);
      for (int k = 0; k < parts; k++)
	derivativeAtQPs[iq] += grdLambda[k][comp] * grd1[k];
    }
  }
1691
1692
1693
1694
1695

  // Some free functions used in MTL4

  template <typename T>
  inline std::size_t size(const AMDiS::DOFVector<T>& v)
Thomas Witkowski's avatar
Thomas Witkowski committed
1696
  {
1697
1698
    return v.getSize();
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1699

1700
1701
1702
1703
1704
1705
  template <typename T>
  inline std::size_t num_rows(const AMDiS::DOFVector<T>& v)
  {
    return v.getSize();
  }

1706

1707
1708
1709
1710
1711
1712
  template <typename T>
  inline std::size_t num_cols(const AMDiS::DOFVector<T>& v)
  {
    return 1;
  }
 
1713

1714
1715
1716
1717
1718
1719
1720
  template <typename T>
  inline void set_to_zero(AMDiS::DOFVector<T>& v)
  {
    using math::zero;
    T  ref, my_zero(zero(ref));

    std::fill(v.begin(), v.end(), my_zero);
Thomas Witkowski's avatar
Thomas Witkowski committed
1721
1722
  }

1723

1724
1725
1726
  template<typename T>
  double DOFVector<T>::DoubleWell(Quadrature* q) const
  {
1727
1728
    Mesh* mesh = this->feSpace->getMesh();

1729
    if (!q) {
1730
      int deg = 2 * this->feSpace->getBasisFcts()->getDegree();
Thomas Witkowski's avatar
Thomas Witkowski committed
1731
      q = Quadrature::provideQuadrature(this->dim, deg);
1732
1733
    }

1734
1735
1736
1737
1738
    FastQuadrature *quadFast = 
      FastQuadrature::provideFastQuadrature(this->feSpace->getBasisFcts(), *q, INIT_PHI);

    double result = 0.0;
    int nPoints = quadFast->getNumPoints();
1739
    mtl::dense_vector<T> uh_vec(nPoints);
1740
1741
    TraverseStack stack;
    ElInfo *elInfo = 
1742
      stack.traverseFirst(mesh, -1,
1743
			  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET);
1744
1745
1746
    while (elInfo) {
      double det = elInfo->getDet();
      double normT = 0.0;
1747
      this->getVecAtQPs(elInfo, nullptr, quadFast, uh_vec);
1748
1749
1750
      for (int iq = 0; iq < nPoints; iq++)
	normT += quadFast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);
      result += det * normT;
1751

1752
1753
      elInfo = stack.traverseNext(elInfo);
    }
1754

1755
1756
1757
1758
1759
1760
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    double localResult = result;
    MPI::COMM_WORLD.Allreduce(&localResult, &result, 1, MPI_DOUBLE, MPI_SUM);
#endif
    
    return result;
1761
  }
1762
1763
1764
1765
1766
1767


  template<typename T>
  DOFVector<typename GradientType<T>::type>*
  DOFVector<T>::getGradient(DOFVector<typename GradientType<T>::type> *grad) const
  {
1768
    FUNCNAME_DBG("DOFVector<T>::getGradient()");
1769
1770
1771
    const FiniteElemSpace *feSpace = DOFVector<T>::feSpace;

    // define result vector
1772
    static DOFVector<typename GradientType<T>::type> *result = nullptr; // TODO: REMOVE STATIC
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825

    if (grad) {
      result = grad;
    } else {
      if(result && result->getFeSpace() != feSpace) {
	delete result;
	result = new DOFVector<typename GradientType<T>::type>(feSpace, "gradient");
      }
    }

    Mesh *mesh = feSpace->getMesh();
    int dim = mesh->getDim();
    const BasisFunction *basFcts = feSpace->getBasisFcts();
    DOFAdmin *admin = feSpace->getAdmin();

    // count number of nodes and dofs per node
    std::vector<int> nNodeDOFs;
    std::vector<int> nNodePreDofs;
    std::vector<DimVec<double>*> bary;

    int nNodes = 0;
    int nDofs = 0;

    for (int i = 0; i < dim + 1; i++) {
      GeoIndex geoIndex = INDEX_OF_DIM(i, dim);
      int nPositions = mesh->getGeo(geoIndex);
      int numPreDofs = admin->getNumberOfPreDofs(i);
      for (int j = 0; j < nPositions; j++) {
	int dofs = basFcts->getNumberOfDofs(geoIndex);
	nNodeDOFs.push_back(dofs);
	nDofs += dofs;
	nNodePreDofs.push_back(numPreDofs);
      }
      nNodes += nPositions;
    }

    TEST_EXIT_DBG(nDofs == basFcts->getNumber())
      ("number of dofs != number of basis functions\n");

    for (int i = 0; i < nDofs; i++)
      bary.push_back(basFcts->getCoords(i));

    mtl::dense_vector<T> localUh(basFcts->getNumber());

    // traverse mesh
    std::vector<bool> visited(getUsedSize(), false);
    TraverseStack stack;
    Flag fillFlag = Mesh::CALL_LEAF_EL | Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);

    while (elInfo) {
      const DegreeOfFreedom **dof = elInfo->getElement()->getDof();
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
1826
      this->getLocalVector(elInfo->getElement(), localUh);
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856

      int localDOFNr = 0;
      for (int i = 0; i < nNodes; i++) { // for all nodes
	for (int j = 0; j < nNodeDOFs[i]; j++) { // for all dofs at this node
	  DegreeOfFreedom dofIndex = dof[i][nNodePreDofs[i] + j];
	  if (!visited[dofIndex]) {
	    basFcts->evalGrdUh(*(bary[localDOFNr]), grdLambda,
			       localUh, ((*result)[dofIndex]));

	    visited[dofIndex] = true;
	  }
	  localDOFNr++;
	}
      }

      elInfo = stack.traverseNext(elInfo);
    }

    return result;
  }


  template<typename T>
  DOFVector<typename GradientType<T>::type>*
  DOFVector<T>::getRecoveryGradient(DOFVector<typename GradientType<T>::type> *grad) const
  {
    const FiniteElemSpace *feSpace = DOFVector<T>::feSpace;
    int dim = DOFVector<T>::dim;

    // define result vector
1857
    static DOFVector<typename GradientType<T>::type> *vec = nullptr; // TODO: REMOVE STATIC
1858
1859
1860
1861
1862
1863

    DOFVector<typename GradientType<T>::type> *result = grad;

    if (!result) {
      if (vec && vec->getFeSpace() != feSpace) {
	delete vec;
1864
	vec = nullptr;
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
      }
      if (!vec)
	vec = new DOFVector<typename GradientType<T>::type>(feSpace, "gradient");

      result = vec;
    }

    typename GradientType<T>::type grd;
    nullify(grd);

    result->set(grd);

    DOFVector<double> volume(feSpace, "volume");
    volume.set(0.0);

    const BasisFunction *basFcts = feSpace->getBasisFcts();
    int nBasisFcts = basFcts->getNumber();
    DimVec<double> bary(dim, DEFAULT_VALUE, (1.0 / (dim + 1.0)));

    // traverse mesh
    Mesh *mesh = feSpace->getMesh();
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh, -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_DET |
					 Mesh::FILL_GRD_LAMBDA | Mesh::FILL_COORDS);

    mtl::dense_vector<T> localUh(nBasisFcts);
    std::vector<DegreeOfFreedom> localIndices(nBasisFcts);

    while (elInfo) {
      double det = elInfo->getDet();
      const DimVec<WorldVector<double> > &grdLambda = elInfo->getGrdLambda();
Naumann, Andreas's avatar
Naumann, Andreas committed
1897
      this->getLocalVector(elInfo->getElement(), localUh);
1898
1899
1900
1901
1902
1903
1904
1905
1906
      basFcts->evalGrdUh(bary, grdLambda, localUh, grd);
      basFcts->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);

      for (int i = 0; i < nBasisFcts; i++) {
	(*result)[localIndices[i]] += grd * det;
	volume[localIndices[i]] += det;
      }

      elInfo = stack.traverseNext(elInfo);
Nestler, Michael's avatar
Nestler, Michael committed
1907
    
1908
    }
Nestler, Michael's avatar
Nestler, Michael committed
1909
1910
    
    
1911
1912
1913
1914
1915
1916
    DOFVector<double>::Iterator volIt(&volume, USED_DOFS);
    typename DOFVector<typename GradientType<T>::type>::Iterator grdIt(result, USED_DOFS);

    for (volIt.reset(), grdIt.reset(); !volIt.end(); ++volIt, ++grdIt)
      if (*volIt != 0.0)
	*grdIt *= 1.0 / (*volIt);
Nestler, Michael's avatar
Nestler, Michael committed
1917
      
1918
1919
1920
1921
1922
1923
1924
1925
    return result;
  }


  template<typename T>
  std::vector<DOFVector<double>*> *transform(DOFVector<typename GradientType<T>::type> *vec,
					     std::vector<DOFVector<double>*> *res)
  {
1926
    FUNCNAME_DBG("DOFVector<T>::transform()");
1927
1928
1929

    TEST_EXIT_DBG(vec)("no vector\n");

1930
    static std::vector<DOFVector<double>*> *result = nullptr; // TODO: REMOVE STATIC
1931
1932
1933
1934
1935
1936
1937

    int len = num_rows(GradientType<T>::getValues((*vec)[0]));
    
    if (!res && !result) {
      result = new std::vector<DOFVector<double>*>(len);
      for (int i = 0; i < len; i++)
	(*result)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
1938
1939
    } else if (res->size() == 0 || (*res)[0] == nullptr) {
      res->resize(len, nullptr);
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
      for (int i = 0; i < len; i++)
	(*res)[i] = new DOFVector<double>(vec->getFeSpace(), "transform");
    }

    std::vector<DOFVector<double>*> *r = res ? res : result;

    int vecSize = vec->getSize();
    for (int i = 0; i < vecSize; i++)
      for (int j = 0; j < len; j++)
	(*((*r)[j]))[i] = (GradientType<T>::getValues((*vec)[i]))[j];

    return r;
  }
1953
}
For faster browsing, not all history is shown. View entire blame