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
1187
    TEST_EXIT(dv)("No vector!\n");

    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());
  }



1188
1189
1190
1191
  void MatrixFct_SOT::initElement(const ElInfo* elInfo, 
				  SubAssembler* subAssembler,
				  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1192
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1193
1194
1195
1196
  }
 
  void VecAtQP_SOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
Thomas Witkowski's avatar
Thomas Witkowski committed
1197
				Quadrature* quad) 
1198
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1199
1200
1201
1202
1203
1204
1205
1206
1207
    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);
1208
  }
1209
1210
1211
1212
1213

  void Vec2AtQP_SOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1214
1215
    vecAtQPs1 = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    vecAtQPs2 = getVectorAtQPs(vec2, elInfo, subAssembler, quad);
1216
  }
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
 
  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
1229
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1230
1231
1232
1233
1234
1235
  }

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

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

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

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

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

  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
1305
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
1306
1307
1308
1309
1310
1311
  }

  void VectorFct_FOT::initElement(const ElInfo* elInfo, 
				  SubAssembler* subAssembler,
				  Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1312
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
  }
 
  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
1326
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
1327
1328
  }

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


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

  void Vec2AtQP_ZOT::initElement(const ElInfo* elInfo, 
				 SubAssembler* subAssembler,
				 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
    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);
1364
1365
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1583
1584
  void VecAtQP_ZOT::eval(int nPoints,
			 const double *uhAtQP,
1585
1586
1587
1588
1589
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
1590
1591
1592
1593
1594
1595
1596
1597
    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];
      }      
1598
1599
1600
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1601

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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



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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
1983
  void Vec2AtQP_SOT::eval(int nPoints,
1984
1985
1986
1987
1988
1989
1990
1991
1992
			  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
1993
      for (int iq = 0; iq < nPoints; iq++) {
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