Operator.cc 95.2 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186

    auxFESpaces.push_back(dv->getFESpace());
  }

  Vec2AtQP_SOT::Vec2AtQP_SOT(DOFVectorBase<double> *dv1, DOFVectorBase<double> *dv2, 
			     BinaryAbstractFunction<double, double, double> *af)
    : SecondOrderTerm(af->getDegree()), vec1(dv1), vec2(dv2), f(af)
  {
    setSymmetric(true);

    TEST_EXIT(dv1)("No first vector!\n");
    TEST_EXIT(dv2)("No second vector!\n");

    auxFESpaces.push_back(dv1->getFESpace());
    auxFESpaces.push_back(dv2->getFESpace());
  }

  MatrixGradient_SOT::MatrixGradient_SOT(DOFVectorBase<double> *dv,
					 AbstractFunction<WorldMatrix<double>, WorldVector<double> > *af,
					 AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
					 bool symm) 
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), divFct(divAf), symmetric(symm)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  FctGradient_SOT::FctGradient_SOT(DOFVectorBase<double> *dv,
				   AbstractFunction<double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndGradientAtQP_SOT::VecAndGradientAtQP_SOT(DOFVectorBase<double> *dv,
						 BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecMatrixGradientAtQP_SOT::VecMatrixGradientAtQP_SOT(DOFVectorBase<double> *dv,
						       BinaryAbstractFunction<WorldMatrix<double>, double, WorldVector<double> > *af,
						       AbstractFunction<WorldVector<double>, WorldMatrix<double> > *divAf,
						       bool symm)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), 
      divFct(divAf), symmetric(symm)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecGradCoordsAtQP_SOT::VecGradCoordsAtQP_SOT(DOFVectorBase<double> *dv,
					       TertiaryAbstractFunction<double, double,
					       WorldVector<double>, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAndCoordsAtQP_SOT::VecAndCoordsAtQP_SOT(DOFVectorBase<double> *dv, 
					     BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af)
  {
    setSymmetric(true);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  MatrixGradientAndCoords_SOT::MatrixGradientAndCoords_SOT(DOFVectorBase<double> *dv,
							   BinaryAbstractFunction<WorldMatrix<double>,
							   WorldVector<double>, WorldVector<double> > *af,
							   AbstractFunction<WorldVector<double>,
							   WorldMatrix<double> > *divAf,
							   bool symm) 
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), divFct(divAf), symmetric(symm)
  {
    setSymmetric(symmetric);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  General_SOT::General_SOT(std::vector<DOFVectorBase<double>*> vecs,
			   std::vector<DOFVectorBase<double>*> grads,
			   TertiaryAbstractFunction<WorldMatrix<double>, 
  			     WorldVector<double>,
			     std::vector<double>, 
			     std::vector<WorldVector<double> > > *f,
			   AbstractFunction<WorldVector<double>, 
			     WorldMatrix<double> > *divFct,
			     bool symmetric)
      : SecondOrderTerm(f->getDegree()),
	vecs_(vecs),
	grads_(grads),
	f_(f),
	divFct_(divFct),
	symmetric_(symmetric)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }
  }

  GeneralParametric_SOT::GeneralParametric_SOT(std::vector<DOFVectorBase<double>*> vecs,
					       std::vector<DOFVectorBase<double>*> grads,
					       QuartAbstractFunction<WorldMatrix<double>, 
					         WorldVector<double>,
					         WorldVector<double>,
					         std::vector<double>, 
					         std::vector<WorldVector<double> > > *f,
					       AbstractFunction<WorldVector<double>, 
					         WorldMatrix<double> > *divFct,
					       bool symmetric)
    : SecondOrderTerm(f->getDegree()),
      vecs_(vecs),
      grads_(grads),
      f_(f),
      divFct_(divFct),
      symmetric_(symmetric)
  {
    vecsAtQPs_.resize(vecs_.size());
    gradsAtQPs_.resize(grads_.size());

    for (int i = 0; i < static_cast<int>(vecs.size()); i++) {
      TEST_EXIT(vecs[i])("One vector is NULL!\n");

      auxFESpaces.push_back(vecs[i]->getFESpace());
    }   

    for (int i = 0; i < static_cast<int>(grads.size()); i++) {
      TEST_EXIT(grads[i])("One gradient vector is NULL!\n");

      auxFESpaces.push_back(grads[i]->getFESpace());
    }
  }

  VecAndGradAtQP_SOT::VecAndGradAtQP_SOT(DOFVectorBase<double> *dv,
					 BinaryAbstractFunction<double, double, WorldVector<double> > *af)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af) 
  {
    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }

  VecAtQP_IJ_SOT::VecAtQP_IJ_SOT(DOFVectorBase<double> *dv, 
				 AbstractFunction<double, double> *af,
				 int x_i, int x_j)
    : SecondOrderTerm(af->getDegree()), vec(dv), f(af), xi(x_i), xj(x_j)
  {
    setSymmetric(xi == xj);

    TEST_EXIT(dv)("No vector!\n");

    auxFESpaces.push_back(dv->getFESpace());
  }



1187
1188
1189
1190
  void MatrixFct_SOT::initElement(const ElInfo* elInfo, 
				  SubAssembler* subAssembler,
				  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1191
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1192
1193
1194
1195
  }
 
  void VecAtQP_SOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
Thomas Witkowski's avatar
Thomas Witkowski committed
1196
				Quadrature* quad) 
1197
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1198
1199
1200
1201
1202
1203
1204
1205
1206
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
  }

  void VecAtQP_SOT::initElement(const ElInfo* smallElInfo, 
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature* quad) 
  {
    vecAtQPs = getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
1207
  }
1208
1209
1210
1211
1212

  void Vec2AtQP_SOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1213
1214
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
1215
  }
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
 
  void CoordsAtQP_SOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void MatrixGradient_SOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1228
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1229
1230
1231
1232
1233
1234
  }

  void FctGradient_SOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1235
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1236
1237
1238
1239
1240
1241
  }

  void VecAndGradientAtQP_SOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1242
1243
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1244
1245
  }

1246
1247
1248
1249
  void VecMatrixGradientAtQP_SOT::initElement(const ElInfo* elInfo, 
					      SubAssembler* subAssembler,
					      Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1250
1251
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1252
1253
  }
  
1254
1255
1256
1257
  void VecAndCoordsAtQP_SOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1258
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1259
1260
1261
1262
1263
1264
1265
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
 
  void MatrixGradientAndCoords_SOT::initElement(const ElInfo* elInfo, 
						SubAssembler* subAssembler,
						Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1266
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1267
1268
1269
1270
1271
1272
1273
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);  
  }

  void VecGradCoordsAtQP_SOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1274
1275
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1276
1277
1278
1279
1280
1281
1282
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);  
  }

  void VecAtQP_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1283
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
  }

  void CoordsAtQP_FOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void VecCoordsAtQP_FOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void VectorGradient_FOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1304
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1305
1306
1307
1308
1309
1310
  }

  void VectorFct_FOT::initElement(const ElInfo* elInfo, 
				  SubAssembler* subAssembler,
				  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1311
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
  }
 
  void VecFctAtQP_FOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) 
  {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

  void VecAtQP_ZOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1325
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1326
1327
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1328
1329
1330
1331
1332
1333
1334
1335
1336
  void VecAtQP_ZOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    vecAtQPs = getVectorAtQPs(vec, smallElInfo, largeElInfo, subAssembler, quad);
  }


1337
1338
1339
1340
  void MultVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1341
1342
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
1343
1344
1345
1346
1347
1348
  }

  void Vec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
  }

  void Vec2AtQP_ZOT::initElement(const ElInfo* smallElInfo, 
				 const ElInfo* largeElInfo,
				 SubAssembler* subAssembler,
				 Quadrature *quad)
  {
    TEST_EXIT(vec1->getFESpace() == vec2->getFESpace())
      ("Not yet implemented!\n");

    vecAtQPs1 = getVectorAtQPs(vec1, smallElInfo, largeElInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad);
1363
1364
  }

1365
1366
1367
1368
  void Vec3AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1369
1370
1371
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
    vecAtQPs3 = getVectorAtQPs(vec3, elInfo, subAssembler, quad);
1372
1373
  }

1374
1375
1376
1377
  void VecAndCoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1378
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1379
1380
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
1381
1382
1383
1384
1385
  
  void Vec2AndGradAtQP_ZOT::initElement(const ElInfo* elInfo, 
					SubAssembler* subAssembler,
					Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1386
1387
1388
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec1, elInfo, subAssembler, quad);
1389
1390
  }

1391
1392
1393
1394
  void FctGradientCoords_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1395
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1396
1397
1398
1399
1400
1401
1402
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
 
  void VecGradCoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1403
1404
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1405
1406
1407
1408
1409
1410
1411
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }
 
  void VecAndGradAtQP_ZOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1412
1413
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1414
1415
1416
1417
1418
1419
  }
 
  void VecAndGradVecAtQP_ZOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1420
1421
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vecGrd, elInfo, subAssembler, quad);
1422
1423
1424
1425
1426
1427
  }

  void VecAndGradVec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1428
1429
1430
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    grad1AtQPs = getGradientsAtQPs(vecGrd1, elInfo, subAssembler, quad);
    grad2AtQPs = getGradientsAtQPs(vecGrd2, elInfo, subAssembler, quad);
1431
1432
1433
1434
1435
1436
  }
 
  void FctGradient_ZOT::initElement(const ElInfo* elInfo, 
				    SubAssembler* subAssembler,
				    Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1437
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1438
1439
1440
1441
1442
1443
1444
1445
  }
 
  void CoordsAtQP_ZOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad) {
    coordsAtQPs = subAssembler->getCoordsAtQPs(elInfo, quad);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1446
  void MatrixFct_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
1447
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1448
    for (int iq = 0; iq < nPoints; iq++)
1449
1450
1451
      lalt(Lambda, (*matrixFct)(vecAtQPs[iq]), *(LALt[iq]), symmetric, 1.0);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1452
  void VecAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
1453
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1454
    for (int iq = 0; iq < nPoints; iq++)
1455
1456
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq]));
  }
1457

Thomas Witkowski's avatar
Thomas Witkowski committed
1458
  void Vec2AtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
1459
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1460
    for (int iq = 0; iq < nPoints; iq++)
1461
1462
1463
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs1[iq], vecAtQPs2[iq]));
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1464
  void CoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
1465
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1466
    for (int iq = 0; iq < nPoints; iq++)
1467
1468
1469
      l1lt(Lambda, (*LALt[iq]), (*g)(coordsAtQPs[iq]));
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1470
  void MatrixGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
1471
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1472
    for (int iq = 0; iq < nPoints; iq++) {
1473
      lalt(Lambda, (*f)(gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
1474
    }
1475
1476
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1477
  void FctGradient_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
1478
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1479
    for (int iq = 0; iq < nPoints; iq++)
1480
1481
1482
      l1lt(Lambda, *(LALt[iq]), (*f)(gradAtQPs[iq]));
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1483
  void VecAndGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
1484
1485
				       DimMat<double> **LALt) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1486
    for (int iq = 0; iq < nPoints; iq++)
1487
1488
1489
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1490
  void VecMatrixGradientAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
1491
1492
					  DimMat<double> **LALt) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1493
    for (int iq = 0; iq < nPoints; iq++)
1494
1495
1496
      lalt(Lambda, (*f)(vecAtQPs[iq], gradAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1497
  void VecAndCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt)
1498
    const { const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1499
    for (int iq = 0; iq < nPoints; iq++)
1500
1501
1502
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], coordsAtQPs[iq]));
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1503
  void MatrixGradientAndCoords_SOT::getLALt(const ElInfo *elInfo, int nPoints,
1504
1505
					    DimMat<double> **LALt) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1506
    for (int iq = 0; iq < nPoints; iq++)
1507
1508
1509
      lalt(Lambda, (*f)(gradAtQPs[iq], coordsAtQPs[iq]), (*LALt[iq]), symmetric, 1.0);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1510
  void VecGradCoordsAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints,
1511
1512
				      DimMat<double> **LALt) const {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
1513
    for (int iq = 0; iq < nPoints; iq++)
1514
1515
1516
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]));
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1517
  void VecAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
1518
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
1519
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1520
      if (b)
1521
1522
1523
1524
1525
1526
	lb(Lambda, *b, Lb[iq], (*f)(vecAtQPs[iq]));
      else
	l1(Lambda, Lb[iq], (*f)(vecAtQPs[iq]));
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1527
  void CoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
1528
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
1529
    for (int iq = 0; iq < nPoints; iq++) {
1530
1531
1532
1533
      l1(Lambda, Lb[iq], (*g)(coordsAtQPs[iq]));
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1534
  void VecCoordsAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints,
1535
1536
1537
				VectorOfFixVecs<DimVec<double> >& Lb) const
  {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
1538
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1539
1540
      lb(Lambda, b, Lb[iq], (*g)(coordsAtQPs[iq]));
    }
1541
1542
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1543
  void VectorGradient_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
1544
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
1545
    if (f) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1546
      for (int iq = 0; iq < nPoints; iq++) {
1547
1548
1549
	lb(Lambda, (*f)(gradAtQPs[iq]), Lb[iq], 1.0);
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1550
      for (int iq = 0; iq < nPoints; iq++) {
1551
1552
1553
1554
1555
	lb(Lambda, gradAtQPs[iq], Lb[iq], 1.0);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1556
  void VectorFct_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
1557
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
1558
    for (int iq = 0; iq < nPoints; iq++) {
1559
1560
1561
1562
      lb(Lambda, (*vecFct)(vecAtQPs[iq]), Lb[iq], 1.0);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1563
  void VecFctAtQP_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
1564
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
1565
    for (int iq = 0; iq < nPoints; iq++) {
1566
1567
1568
1569
      lb(Lambda, (*g)(coordsAtQPs[iq]), Lb[iq], 1.0);
    }
  }

1570
  void VecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1571
    if (f) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1572
      for (int iq = 0; iq < nPoints; iq++) {
1573
1574
1575
	C[iq] += (*f)(vecAtQPs[iq]);
      }
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1576
      for (int iq = 0; iq < nPoints; iq++) {
1577
1578
1579
1580
1581
	C[iq] += vecAtQPs[iq];
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1582
1583
  void VecAtQP_ZOT::eval(int nPoints,
			 const double *uhAtQP,
1584
1585
1586
1587
1588
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1589
1590
1591
1592
1593
1594
1595
1596
    if (f) {
      for (int iq = 0; iq < nPoints; iq++) {
	result[iq] += fac * (*f)(vecAtQPs[iq]) * uhAtQP[iq];
      }
    } else {
      for (int iq = 0; iq < nPoints; iq++) {
	result[iq] += fac * vecAtQPs[iq] * uhAtQP[iq];
      }      
1597
1598
1599
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1600

1601
  void MultVecAtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1602
    for (int iq = 0; iq < nPoints; iq++) {
1603
1604
1605
1606
      C[iq] += (*f1)(vecAtQPs1[iq]) * (*f2)(vecAtQPs2[iq]);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1607
1608
  void MultVecAtQP_ZOT::eval(int nPoints,
			     const double *uhAtQP,
1609
1610
1611
1612
1613
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
			     double *result,
			     double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1614
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1615
1616
      result[iq] += fac * (*f1)(vecAtQPs1[iq]) * 
	(*f2)(vecAtQPs2[iq]) * uhAtQP[iq];
1617
1618
1619
    }
  }

1620
  void Vec2AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1621
    for (int iq = 0; iq < nPoints; iq++) {
1622
1623
1624
1625
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1626
1627
  void Vec2AtQP_ZOT::eval(int nPoints,
			  const double *uhAtQP,
1628
1629
1630
1631
1632
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
			  double *result,
			  double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1633
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1634
      result[iq] += fac * (*f)(vecAtQPs1[iq], vecAtQPs2[iq]) * 
1635
1636
1637
1638
	uhAtQP[iq];
    }
  }

1639
  void Vec3AtQP_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1640
    for (int iq = 0; iq < nPoints; iq++) {
1641
1642
1643
1644
      C[iq] += (*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1645
  void Vec3AtQP_ZOT::eval(int nPoints,
1646
1647
1648
1649
1650
1651
			  const double *uhAtQP,
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
			  double *result,
			  double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1652
    for (int iq = 0; iq < nPoints; iq++) {
1653
1654
1655
1656
1657
1658
1659
1660
      result[iq] += 
	fac * 
	(*f)(vecAtQPs1[iq], vecAtQPs2[iq], vecAtQPs3[iq]) * 
	uhAtQP[iq];
    }
  }


1661
1662
  void VecAndCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints, 
				  std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1663
    for (int iq = 0; iq < nPoints; iq++) {
1664
1665
1666
1667
      C[iq] += (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1668
1669
  void VecAndCoordsAtQP_ZOT::eval(int nPoints,
				  const double *uhAtQP,
1670
1671
1672
1673
1674
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
				  double *result,
				  double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1675
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1676
      result[iq] += fac * (*f)(vecAtQPs[iq], coordsAtQPs[iq]) * 
1677
1678
1679
1680
1681
	uhAtQP[iq];
    }
  }

 
1682
1683
  void FctGradientCoords_ZOT::getC(const ElInfo *, int nPoints, 
				   std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1684
    for (int iq = 0; iq < nPoints; iq++) {
1685
1686
1687
1688
      C[iq] += (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
    }
  }
 
Thomas Witkowski's avatar
Thomas Witkowski committed
1689
1690
  void FctGradientCoords_ZOT::eval(int nPoints,
				   const double *uhAtQP,
1691
1692
1693
1694
1695
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1696
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1697
      result[iq] += fac * (*f)(gradAtQPs[iq], coordsAtQPs[iq]) * 
1698
1699
1700
1701
	uhAtQP[iq];
    }
  }

1702
1703
  void VecGradCoordsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				   std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1704
    for (int iq = 0; iq < nPoints; iq++) {
1705
1706
1707
1708
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
    }
  }
 
Thomas Witkowski's avatar
Thomas Witkowski committed
1709
1710
  void VecGradCoordsAtQP_ZOT::eval(int nPoints,
				   const double *uhAtQP,
1711
1712
1713
1714
1715
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
1716
    for (int iq = 0; iq < nPoints; iq++)
1717
1718
1719
1720
1721
1722
      result[iq] += 
	fac * 
	(*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]) * 
	uhAtQP[iq];
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1723
  void Vec2AndGradAtQP_ZOT::eval(int nPoints,
1724
1725
1726
1727
1728
1729
				 const double *uhAtQP,
				 const WorldVector<double> *grdUhAtQP,
				 const WorldMatrix<double> *D2UhAtQP,
				 double *result,
				 double fac) const
  {
1730
    for (int iq = 0; iq < nPoints; iq++)
1731
1732
1733
1734
1735
1736
      result[iq] += 
	fac * 
	(*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]) * 
	uhAtQP[iq];
  }

1737
1738
  void Vec2AndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
				 std::vector<double> &C) const 
1739
  { 
1740
    for (int iq = 0; iq < nPoints; iq++)
1741
1742
1743
      C[iq] += (*f)(vecAtQPs1[iq], gradAtQPs[iq], vecAtQPs2[iq]);
  }

1744
  void VecAndGradAtQP_ZOT::getC(const ElInfo *, int nPoints,
1745
1746
1747
				std::vector<double> &C) const 
  { 
    for (int iq = 0; iq < nPoints; iq++)
1748
1749
1750
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
  }
 
Thomas Witkowski's avatar
Thomas Witkowski committed
1751
1752
  void VecAndGradAtQP_ZOT::eval(int nPoints,
				const double *uhAtQP,
1753
1754
1755
1756
1757
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1758
    for (int iq = 0; iq < nPoints; iq++)
1759
      result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];   
1760
1761
  }

1762
1763
  void VecAndGradVecAtQP_ZOT::getC(const ElInfo *, int nPoints,
				   std::vector<double> &C) const 
1764
  { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1765
    for (int iq = 0; iq < nPoints; iq++)
1766
1767
1768
      C[iq] += (*f)(vecAtQPs[iq], gradAtQPs[iq]);
  }
 
Thomas Witkowski's avatar
Thomas Witkowski committed
1769
1770
  void VecAndGradVecAtQP_ZOT::eval(int nPoints,
				   const double *uhAtQP,
1771
1772
1773
1774
1775
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1776
1777
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(vecAtQPs[iq], gradAtQPs[iq]) * uhAtQP[iq];
1778
1779
  }

1780
1781
  void VecAndGradVec2AtQP_ZOT::getC(const ElInfo *, int nPoints,
				    std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1782
    for (int iq = 0; iq < nPoints; iq++)
1783
1784
1785
      C[iq] += (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]);
  }
 
Thomas Witkowski's avatar
Thomas Witkowski committed
1786
1787
  void VecAndGradVec2AtQP_ZOT::eval(int nPoints,
				    const double *uhAtQP,
1788
1789
1790
1791
1792
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
				    double *result,
				    double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1793
    for (int iq = 0; iq < nPoints; iq++)
1794
      result[iq] += 
Thomas Witkowski's avatar
Thomas Witkowski committed
1795
	fac * (*f)(vecAtQPs[iq], grad1AtQPs[iq], grad2AtQPs[iq]) * uhAtQP[iq];  
1796
1797
  }

1798
  void FctGradient_ZOT::getC(const ElInfo *, int nPoints, std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1799
    for (int iq = 0; iq < nPoints; iq++)
1800
1801
1802
      C[iq] += (*f)(gradAtQPs[iq]);
  }
 
Thomas Witkowski's avatar
Thomas Witkowski committed
1803
1804
  void FctGradient_ZOT::eval(int nPoints,
			     const double *uhAtQP,
1805
1806
1807
1808
1809
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
			     double *result,
			     double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1810
1811
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*f)(gradAtQPs[iq]) * uhAtQP[iq];    
1812
1813
  }

1814
1815
  void CoordsAtQP_ZOT::getC(const ElInfo *elInfo, int nPoints, 
			    std::vector<double> &C) const { 
Thomas Witkowski's avatar
Thomas Witkowski committed
1816
1817
    for (int iq = 0; iq < nPoints; iq++)
      C[iq] += (*g)(coordsAtQPs[iq]);    
1818
1819
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1820
1821
  void CoordsAtQP_ZOT::eval(int nPoints,
			    const double *uhAtQP,
1822
1823
1824
1825
1826
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1827
1828
    for (int iq = 0; iq < nPoints; iq++)
      result[iq] += fac * (*g)(coordsAtQPs[iq]) * uhAtQP[iq];    
1829
1830
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1831
1832
  void MatrixFct_SOT::eval(int nPoints,
			   const double *uhAtQP,
1833
1834
1835
1836
1837
			   const WorldVector<double> *grdUhAtQP,
			   const WorldMatrix<double> *D2UhAtQP,
			   double *result,
			   double factor) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1838
    int dow = Global::getGeo(WORLD);
1839

Thomas Witkowski's avatar
Thomas Witkowski committed
1840
    for (int iq = 0; iq < nPoints; iq++) {
1841
1842
1843
1844
      double resultQP = 0.0;

      WorldMatrix<double> A = (*matrixFct)(vecAtQPs[iq]);

Thomas Witkowski's avatar
Thomas Witkowski committed
1845
1846
1847
      if (D2UhAtQP)
	for (int i = 0; i < dow; i++)
	  for (int j = 0; j < dow; j++)
1848
1849
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];

Thomas Witkowski's avatar
Thomas Witkowski committed
1850
      if (grdUhAtQP)
1851
1852
1853
1854
1855
1856
	resultQP += (*divFct)(A) * grdUhAtQP[iq];

      result[iq] += resultQP * factor;
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1857
  void MatrixFct_SOT::weakEval(int nPoints,
1858
1859
1860
			       const WorldVector<double> *grdUhAtQP,
			       WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1861
    if (grdUhAtQP) {
1862
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
1863
      for (int iq = 0; iq < nPoints; iq++) {
1864
1865
1866
1867
1868
1869
1870
	A = (*matrixFct)(vecAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
1871
1872
  void Matrix_SOT::eval(int nPoints,
			const double *uhAtQP,
1873
1874
1875
1876
1877
			const WorldVector<double> *grdUhAtQP,
			const WorldMatrix<double> *D2UhAtQP,
			double *result,
			double factor) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1878
    int dow = Global::getGeo(WORLD);
1879
  
Thomas Witkowski's avatar
Thomas Witkowski committed
1880
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1881
      for (int iq = 0; iq < nPoints; iq++) {
1882
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1883
1884
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
1885
1886
1887
1888
1889
1890
1891
1892
	    resultQP += matrix[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
	result[iq] += resultQP * factor;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1893
  void Matrix_SOT::weakEval(int nPoints,
1894
1895
1896
			    const WorldVector<double> *grdUhAtQP,
			    WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1897
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1898
      for (int iq = 0; iq < nPoints; iq++) {
1899
1900
1901
1902
1903
1904
1905
	result[iq] += matrix * grdUhAtQP[iq];
      }
    }
  }



Thomas Witkowski's avatar
Thomas Witkowski committed
1906
1907
  void MatrixGradient_SOT::eval(int nPoints,
				const double *uhAtQP,
1908
1909
1910
1911
1912
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double factor) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1913
    int dow = Global::getGeo(WORLD);
1914

Thomas Witkowski's avatar
Thomas Witkowski committed
1915
    for (int iq = 0; iq < nPoints; iq++) {
1916
1917
1918
1919
      double resultQP = 0.0;

      WorldMatrix<double> A = (*f)(gradAtQPs[iq]);

Thomas Witkowski's avatar
Thomas Witkowski committed
1920
1921
1922
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
1923
1924
1925
1926
1927
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
1928
      if (grdUhAtQP) {
1929
1930
1931
1932
1933
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }

      result[iq] += resultQP * factor;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1934
  }
1935

Thomas Witkowski's avatar
Thomas Witkowski committed
1936
  void MatrixGradient_SOT::weakEval(int nPoints,
1937
1938
1939
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1940
    if (grdUhAtQP) {
1941
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
1942
      for (int iq = 0; iq < nPoints; iq++) {
1943
1944
1945
1946
1947
1948
	A = (*f)(gradAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1949
1950
  void VecAtQP_SOT::eval(int nPoints,
			 const double *uhAtQP,
1951
1952
1953
1954
1955
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1956
    int dow = Global::getGeo(WORLD);
1957

Thomas Witkowski's avatar
Thomas Witkowski committed
1958
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1959
      for (int iq = 0; iq < nPoints; iq++) {
1960
1961
	double factor = (*f)(vecAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1962
	for (int i = 0; i < dow; i++) {
1963
1964
1965
1966
1967
1968
1969
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1970
  void VecAtQP_SOT::weakEval(int nPoints,
1971
1972
1973
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1974
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1975
      for (int iq = 0; iq < nPoints; iq++) {
1976
1977
1978
1979
1980
1981
	double factor = (*f)(vecAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1982
  void Vec2AtQP_SOT::eval(int nPoints,
1983
1984
1985
1986
1987
1988
1989
1990
1991
			  const double *uhAtQP,
			  const WorldVector<double> *grdUhAtQP,
			  const WorldMatrix<double> *D2UhAtQP,
			  double *result,
			  double fac) const
  {
    int dow = Global::getGeo(WORLD);
    
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1992
      for (int iq = 0; iq < nPoints; iq++) {
1993
1994
1995
1996
1997
1998
1999
2000
	double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++) {
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
For faster browsing, not all history is shown. View entire blame