Operator.cc 95.2 KB
Newer Older
2001
2002
2003
    }
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
2004
  void Vec2AtQP_SOT::weakEval(int nPoints,
2005
2006
2007
2008
			      const WorldVector<double> *grdUhAtQP,
			      WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2009
      for (int iq = 0; iq < nPoints; iq++) {
2010
2011
2012
2013
2014
2015
	double factor = (*f)(vecAtQPs1[iq], vecAtQPs2[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2016
2017
  void CoordsAtQP_SOT::eval(int nPoints,
			    const double *uhAtQP,
2018
2019
2020
2021
2022
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double f) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2023
    int dow = Global::getGeo(WORLD);
2024

Thomas Witkowski's avatar
Thomas Witkowski committed
2025
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2026
      for (int iq = 0; iq < nPoints; iq++) {
2027
2028
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2029
	for (int i = 0; i < dow; i++) {
2030
2031
2032
2033
2034
2035
2036
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += factor * f * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2037
  void CoordsAtQP_SOT::weakEval(int nPoints,
2038
2039
2040
				const WorldVector<double> *grdUhAtQP,
				WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2041
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2042
      for (int iq = 0; iq < nPoints; iq++) {
2043
2044
2045
2046
2047
2048
	double factor = (*g)(coordsAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2049
2050
  void FctGradient_SOT::eval(int nPoints,
			     const double *uhAtQP,
2051
2052
2053
2054
2055
			     const WorldVector<double> *grdUhAtQP,
			     const WorldMatrix<double> *D2UhAtQP,
			     double *result,
			     double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2056
    int dow = Global::getGeo(WORLD);
2057

Thomas Witkowski's avatar
Thomas Witkowski committed
2058
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2059
      for (int iq = 0; iq < nPoints; iq++) {
2060
2061
	double factor = (*f)(gradAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2062
	for (int i = 0; i < dow; i++) {
2063
2064
2065
2066
2067
2068
2069
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2070
  void FctGradient_SOT::weakEval(int nPoints,
2071
2072
2073
				 const WorldVector<double> *grdUhAtQP,
				 WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2074
    if (grdUhAtQP) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
2075
      for (int iq = 0; iq < nPoints; iq++) {
2076
2077
2078
2079
2080
2081
	double factor = (*f)(gradAtQPs[iq]);
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2082
2083
  void VecAndGradientAtQP_SOT::eval(int nPoints,
				    const double *uhAtQP,
2084
2085
2086
2087
2088
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
				    double *result,
				    double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2089
    int dow = Global::getGeo(WORLD);
2090

Thomas Witkowski's avatar
Thomas Witkowski committed
2091
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2092
      for (int iq = 0; iq < nPoints; iq++) {
2093
2094
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2095
	for (int i = 0; i < dow; i++) {
2096
2097
2098
2099
2100
2101
2102
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2103
  void VecAndGradientAtQP_SOT::weakEval(int nPoints,
2104
2105
2106
					const WorldVector<double> *grdUhAtQP,
					WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2107
    if (grdUhAtQP) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
2108
      for (int iq = 0; iq < nPoints; iq++) {
2109
2110
2111
2112
2113
2114
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2115
  void VecMatrixGradientAtQP_SOT::eval(int nPoints,
2116
2117
2118
2119
2120
2121
2122
2123
				       const double *uhAtQP,
				       const WorldVector<double> *grdUhAtQP,
				       const WorldMatrix<double> *D2UhAtQP,
				       double *result,
				       double factor) const
  {
    int dow = Global::getGeo(WORLD);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
2124
    for (int iq = 0; iq < nPoints; iq++) {
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
      double resultQP = 0.0;
      
      WorldMatrix<double> A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
      
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
      }
      
      if (grdUhAtQP) {
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }
      
      result[iq] += resultQP * factor;
    }
  }
2144

Thomas Witkowski's avatar
Thomas Witkowski committed
2145
  void VecMatrixGradientAtQP_SOT::weakEval(int nPoints,
2146
2147
2148
2149
2150
					   const WorldVector<double> *grdUhAtQP,
					   WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
2151
      for (int iq = 0; iq < nPoints; iq++) {
2152
2153
2154
2155
2156
	A = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }
2157

Thomas Witkowski's avatar
Thomas Witkowski committed
2158
2159
  void VecAndCoordsAtQP_SOT::eval(int nPoints,
				  const double *uhAtQP,
2160
2161
2162
2163
2164
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
				  double *result,
				  double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2165
    int dow = Global::getGeo(WORLD);
2166

Thomas Witkowski's avatar
Thomas Witkowski committed
2167
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2168
      for (int iq = 0; iq < nPoints; iq++) {
2169
2170
	double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2171
	for (int i = 0; i < dow; i++) {
2172
2173
2174
2175
2176
2177
2178
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2179
  void VecAndCoordsAtQP_SOT::weakEval(int nPoints,
2180
2181
2182
				      const WorldVector<double> *grdUhAtQP,
				      WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2183
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2184
      for (int iq = 0; iq < nPoints; iq++) {
2185
2186
2187
2188
2189
2190
	double factor = (*f)(vecAtQPs[iq], coordsAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2191
2192
  void MatrixGradientAndCoords_SOT::eval(int nPoints,
					 const double *uhAtQP,
2193
2194
2195
2196
2197
					 const WorldVector<double> *grdUhAtQP,
					 const WorldMatrix<double> *D2UhAtQP,
					 double *result,
					 double factor) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2198
    int dow = Global::getGeo(WORLD);
2199

Thomas Witkowski's avatar
Thomas Witkowski committed
2200
    for (int iq = 0; iq < nPoints; iq++) {
2201
2202
2203
2204
      double resultQP = 0.0;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2205
2206
2207
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
      }

      if(grdUhAtQP) {
	resultQP += (*divFct)(A) * grdUhAtQP[iq];
      }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2221
  void MatrixGradientAndCoords_SOT::weakEval(int nPoints,
2222
2223
2224
					     const WorldVector<double> *grdUhAtQP,
					     WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2225
    if (grdUhAtQP) {
2226
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
2227
      for (int iq = 0; iq < nPoints; iq++) {
2228
2229
2230
2231
2232
2233
	A = (*f)(gradAtQPs[iq], coordsAtQPs[iq]);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2234
2235
  void VecGradCoordsAtQP_SOT::eval(int nPoints,
				   const double *uhAtQP,
2236
2237
2238
2239
2240
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2241
    int  dow = Global::getGeo(WORLD);
2242

Thomas Witkowski's avatar
Thomas Witkowski committed
2243
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2244
      for (int iq = 0; iq < nPoints; iq++) {
2245
2246
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2247
	for (int i = 0; i < dow; i++) {
2248
2249
2250
2251
2252
2253
2254
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2255
  void VecGradCoordsAtQP_SOT::weakEval(int nPoints,
2256
2257
2258
				       const WorldVector<double> *grdUhAtQP,
				       WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2259
    if (grdUhAtQP) {  
Thomas Witkowski's avatar
Thomas Witkowski committed
2260
      for (int iq = 0; iq < nPoints; iq++) {
2261
2262
2263
2264
2265
2266
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq], coordsAtQPs[iq]);
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2267
2268
  void VecAtQP_FOT::eval(int nPoints,
			 const double *uhAtQP,
2269
2270
2271
2272
2273
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2274
    int dow = Global::getGeo(WORLD);
2275
  
Thomas Witkowski's avatar
Thomas Witkowski committed
2276
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2277
      for (int iq = 0; iq < nPoints; iq++) {
2278
2279
	double factor = (*f)(vecAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2280
	for (int i = 0; i < dow; i++) {
2281
2282
2283
2284
2285
2286
2287
	  resultQP += grdUhAtQP[iq][i];
	}
	result[iq] += fac * factor * resultQP;
      } 
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2288
2289
  void CoordsAtQP_FOT::eval(int nPoints,
			    const double *uhAtQP,
2290
2291
2292
2293
2294
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double f) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2295
    int dow = Global::getGeo(WORLD);
2296
  
Thomas Witkowski's avatar
Thomas Witkowski committed
2297
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2298
      for (int iq = 0; iq < nPoints; iq++) {
2299
2300
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2301
	for (int i = 0; i < dow; i++) {
2302
2303
2304
2305
2306
2307
2308
	  resultQP += grdUhAtQP[iq][i];
	}
	result[iq] += f * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2309
2310
  void VecCoordsAtQP_FOT::eval(int nPoints,
			       const double *uhAtQP,
2311
2312
2313
2314
2315
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
			       double *result,
			       double f) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2316
    int dow = Global::getGeo(WORLD);
2317
  
Thomas Witkowski's avatar
Thomas Witkowski committed
2318
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2319
      for (int iq = 0; iq < nPoints; iq++) {
2320
2321
	double factor = (*g)(coordsAtQPs[iq]);
	double resultQP = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
2322
	for (int i = 0; i < dow; i++) {
2323
2324
2325
2326
2327
2328
2329
	  resultQP += grdUhAtQP[iq][i];
	}
	result[iq] += f * factor * resultQP;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2330
2331
  void VectorGradient_FOT::eval(int nPoints,
				const double *uhAtQP,
2332
2333
2334
2335
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double factor) const
Thomas Witkowski's avatar
Thomas Witkowski committed
2336
2337
2338
  { 
    if (grdUhAtQP) {
      if (f) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2339
	for (int iq = 0; iq < nPoints; iq++) {
2340
2341
2342
2343
	  WorldVector<double> b = (*f)(gradAtQPs[iq]);
	  result[iq] += b * grdUhAtQP[iq] * factor;
	}
      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
2344
	for (int iq = 0; iq < nPoints; iq++) {
2345
2346
2347
2348
2349
2350
	  result[iq] += gradAtQPs[iq] * grdUhAtQP[iq] * factor;
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2351
2352
  void VectorFct_FOT::eval(int nPoints,
			   const double *uhAtQP,
2353
2354
2355
2356
2357
			   const WorldVector<double> *grdUhAtQP,
			   const WorldMatrix<double> *D2UhAtQP,
			   double *result,
			   double factor) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2358
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2359
      for (int iq = 0; iq < nPoints; iq++) {
2360
2361
2362
2363
2364
2365
	WorldVector<double> b = (*vecFct)(vecAtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2366
2367
  void VecFctAtQP_FOT::eval(int nPoints,
			    const double *uhAtQP,
2368
2369
2370
2371
2372
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double f) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2373
    int dow = Global::getGeo(WORLD);
2374
  
Thomas Witkowski's avatar
Thomas Witkowski committed
2375
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2376
      for (int iq = 0; iq < nPoints; iq++) {
2377
2378
	double resultQP = 0.0;
	const WorldVector<double> &b = (*g)(coordsAtQPs[iq]);
Thomas Witkowski's avatar
Thomas Witkowski committed
2379
	for (int i = 0; i < dow; i++) {
2380
2381
2382
2383
2384
2385
2386
	  resultQP += b[i] * grdUhAtQP[iq][i];
	}
	result[iq] += f * resultQP;
      }
    }
  }

2387
  Assembler *Operator::getAssembler(int rank) 
2388
  {
2389
2390
    if (!assembler[rank]) {
      initAssembler(rank, NULL, NULL, NULL, NULL);
2391
2392
    }

2393
2394
2395
2396
2397
2398
    return assembler[rank];
  }

  void Operator::setAssembler(int rank, Assembler *a)
  {
    assembler[rank] = a; 
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
  }

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

  void CoordsAtQP_IJ_SOT::getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
2409
				  int nPoints, 
2410
2411
2412
				  DimMat<double> **LALt) const
  {
    const DimVec<WorldVector<double> > Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
2413
    for (int iq = 0; iq < nPoints; iq++) {
2414
2415
2416
2417
      lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*g)(coordsAtQPs[iq]));
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2418
2419
  void CoordsAtQP_IJ_SOT::eval(int nPoints,
			       const double *uhAtQP,
2420
2421
2422
2423
2424
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
			       double *result,
			       double f) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2425
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2426
      for (int iq = 0; iq < nPoints; iq++) {
2427
2428
2429
2430
2431
2432
	double factor = (*g)(coordsAtQPs[iq]);
	result[iq] = D2UhAtQP[iq][xi][xj] * factor * f;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2433
  void CoordsAtQP_IJ_SOT::weakEval(int nPoints,
2434
2435
2436
				   const WorldVector<double> *grdUhAtQP,
				   WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2437
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2438
      for (int iq = 0; iq < nPoints; iq++) {
2439
2440
2441
2442
2443
2444
2445
2446
2447
2448
	double factor = (*g)(coordsAtQPs[iq]);  
	result[iq][xi] += grdUhAtQP[iq][xj] * factor;
      }
    }
  }

  void VecAtQP_IJ_SOT::initElement(const ElInfo* elInfo, 
				   SubAssembler* subAssembler,
				   Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2449
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
2450
2451
2452
  }
 
  void VecAtQP_IJ_SOT::getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
2453
			       int nPoints, 
2454
2455
2456
			       DimMat<double> **LALt) const
  {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
2457
    for (int iq = 0; iq < nPoints; iq++) {
2458
2459
2460
2461
      lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*f)(vecAtQPs[iq]));
    }
  }  

Thomas Witkowski's avatar
Thomas Witkowski committed
2462
2463
  void VecAtQP_IJ_SOT::eval(int nPoints,
			    const double *uhAtQP,
2464
2465
2466
2467
2468
			    const WorldVector<double> *grdUhAtQP,
			    const WorldMatrix<double> *D2UhAtQP,
			    double *result,
			    double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2469
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2470
      for (int iq = 0; iq < nPoints; iq++) {
2471
2472
2473
2474
2475
2476
	double factor = (*f)(vecAtQPs[iq]);
	result[iq] = D2UhAtQP[iq][xi][xj] * factor * fac;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2477
  void VecAtQP_IJ_SOT::weakEval(int nPoints,
2478
2479
2480
				const WorldVector<double> *grdUhAtQP,
				WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2481
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2482
      for (int iq = 0; iq < nPoints; iq++) {
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
2493
	double factor = (*f)(vecAtQPs[iq]);  
	result[iq][xi] += grdUhAtQP[iq][xj] * factor;
      }
    }
  }

 
  void VecOfDOFVecsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					 SubAssembler* subAssembler,
					 Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2494
2495
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2496
      vecsAtQPs[i] = getVectorAtQPs(vecs[i], elInfo, subAssembler, quad);
2497
2498
2499
    }
  }
 
2500
2501
  void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				  std::vector<double> &C) const 
2502
  { 
2503
    int size = static_cast<int>(vecs.size());
2504
    std::vector<double> arg(size);
2505

Thomas Witkowski's avatar
Thomas Witkowski committed
2506
    for (int iq = 0; iq < nPoints; iq++) {
2507
      for (int i = 0; i < size; i++) {
2508
2509
2510
2511
2512
2513
	arg[i] = vecsAtQPs[i][iq];
      }
      C[iq] += (*f)(arg);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2514
2515
  void VecOfDOFVecsAtQP_ZOT::eval(int nPoints,
				  const double *uhAtQP,
2516
2517
2518
2519
2520
				  const WorldVector<double> *grdUhAtQP,
				  const WorldMatrix<double> *D2UhAtQP,
				  double *result,
				  double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2521
    int size = static_cast<int>(vecs.size());
2522
    std::vector<double> arg(size);
2523

Thomas Witkowski's avatar
Thomas Witkowski committed
2524
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2525
      for (int i = 0; i < size; i++) {
2526
2527
2528
2529
2530
2531
2532
2533
2534
2535
2536
2537
2538
	arg[i] = vecsAtQPs[i][iq];
      }
      result[iq] += fac * (*f)(arg) * uhAtQP[iq];
    }
  }

 
  void VecOfGradientsAtQP_ZOT::initElement(const ElInfo* elInfo, 
					   SubAssembler* subAssembler,
					   Quadrature *quad) 
  {
    int i ,size = static_cast<int>(vecs.size());
    for(i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2539
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
2540
2541
2542
2543
2544
2545
2546
    }
  }
 
  void VecDivergence_ZOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2547
2548
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2549
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
2550
2551
2552
2553
    }
  }
 

2554
2555
  void VecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				    std::vector<double> &C) const 
2556
  { 
2557
    int size = static_cast<int>(vecs.size());
2558
    std::vector<WorldVector<double>*> arg(size);
2559

Thomas Witkowski's avatar
Thomas Witkowski committed
2560
    for (int iq = 0; iq < nPoints; iq++) {
2561
      for (int i = 0; i < size; i++) {
2562
2563
2564
2565
2566
2567
	arg[i] = &(gradsAtQPs[i][iq]);
      }
      C[iq] += (*f)(arg);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2568
2569
  void VecOfGradientsAtQP_ZOT::eval(int nPoints,
				    const double *uhAtQP,
2570
2571
2572
2573
2574
				    const WorldVector<double> *grdUhAtQP,
				    const WorldMatrix<double> *D2UhAtQP,
				    double *result,
				    double fac) const
  {
2575
    int size = static_cast<int>(vecs.size());
2576
    std::vector<WorldVector<double>*> arg(size);
2577

Thomas Witkowski's avatar
Thomas Witkowski committed
2578
    for (int iq = 0; iq < nPoints; iq++) {
2579
      for (int i = 0; i < size; i++) {
2580
2581
2582
2583
2584
2585
	arg[i] = &(gradsAtQPs[i][iq]);
      }
      result[iq] += fac * (*f)(arg) * uhAtQP[iq];
    }
  }

2586
2587
  void VecDivergence_ZOT::getC(const ElInfo *, int nPoints,
			       std::vector<double> &C) const 
2588
  { 
2589
    int size = static_cast<int>(vecs.size());
2590

Thomas Witkowski's avatar
Thomas Witkowski committed
2591
    for (int iq = 0; iq < nPoints; iq++) {
2592
      for (int i = 0; i < size; i++) {
2593
2594
2595
2596
2597
	C[iq] += gradsAtQPs[i][iq][i];
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2598
2599
  void VecDivergence_ZOT::eval(int nPoints,
			       const double *uhAtQP,
2600
2601
2602
2603
2604
			       const WorldVector<double> *grdUhAtQP,
			       const WorldMatrix<double> *D2UhAtQP,
			       double *result,
			       double fac) const
  {
2605
    int size = static_cast<int>(vecs.size());
2606

Thomas Witkowski's avatar
Thomas Witkowski committed
2607
    for (int iq = 0; iq < nPoints; iq++) {
2608
2609
      double d = 0.0;
      for (int i = 0; i < size; i++) {
2610
2611
2612
2613
2614
2615
	d += gradsAtQPs[i][iq][i];
      }
      result[iq] += d * uhAtQP[iq] * fac;
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2616
2617
  void VecAndGradAtQP_SOT::eval(int nPoints,
				const double *uhAtQP,
2618
2619
2620
2621
2622
				const WorldVector<double> *grdUhAtQP,
				const WorldMatrix<double> *D2UhAtQP,
				double *result,
				double fac) const
  {
2623
    int dow = Global::getGeo(WORLD);
2624

2625
    if (D2UhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2626
      for (int iq = 0; iq < nPoints; iq++) {
2627
2628
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	double resultQP = 0.0;
2629
	for (int i = 0; i < dow; i++) {
2630
2631
2632
2633
2634
2635
2636
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * resultQP * factor;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2637
  void VecAndGradAtQP_SOT::weakEval(int nPoints,
2638
2639
2640
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
2641
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2642
      for (int iq = 0; iq < nPoints; iq++) {
2643
2644
2645
2646
2647
2648
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2649
  void VecAndGradAtQP_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
2650
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
2651
    for (int iq = 0; iq < nPoints; iq++) {
2652
2653
2654
2655
2656
2657
2658
2659
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
    }
  }
  
  void VecAndGradAtQP_SOT::initElement(const ElInfo* elInfo, 
				       SubAssembler* subAssembler,
				       Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2660
2661
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec, elInfo, subAssembler, quad);
2662
2663
2664
2665
2666
2667
  }
 
  void General_SOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2668
2669
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2670
2671
2672

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

Thomas Witkowski's avatar
Thomas Witkowski committed
2673
2674
    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
2675
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
2676
2677
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
2678
2679
2680
2681
    }
  }

  void General_SOT::getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
2682
			    int nPoints, 
2683
2684
			    DimMat<double> **LALt) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2685
2686
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2687

Thomas Witkowski's avatar
Thomas Witkowski committed
2688
2689
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2690
2691

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
2692

Thomas Witkowski's avatar
Thomas Witkowski committed
2693
2694
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
2695
2696
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2697
      for (int i = 0; i < nGrads; i++) {
2698
2699
2700
2701
2702
2703
2704
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lalt(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), 
	   *(LALt[iq]), symmetric_, 1.0);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2705
2706
  void General_SOT::eval(int nPoints,
			 const double *uhAtQP,
2707
2708
2709
2710
2711
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
2712
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
2713
2714
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2715

Thomas Witkowski's avatar
Thomas Witkowski committed
2716
2717
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2718

Thomas Witkowski's avatar
Thomas Witkowski committed
2719
    for (int iq = 0; iq < nPoints; iq++) {
2720
2721
      double resultQP = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
2722
      for (int i = 0; i < nVecs; i++) {
2723
2724
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2725
      for (int i = 0; i < nGrads; i++) {
2726
2727
2728
2729
2730
	gradsArg[i] = gradsAtQPs_[i][iq];
      }

      WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);

2731
2732
2733
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
2734
2735
2736
2737
2738
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
      }

2739
      if (grdUhAtQP) {
2740
2741
2742
2743
2744
2745
2746
	resultQP += (*divFct_)(A) * grdUhAtQP[iq];
      }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2747
  void General_SOT::weakEval(int nPoints,
2748
2749
2750
			     const WorldVector<double> *grdUhAtQP,
			     WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2751
2752
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2753

Thomas Witkowski's avatar
Thomas Witkowski committed
2754
2755
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2756

2757
    if (grdUhAtQP) {
2758
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
2759
2760
      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nVecs; i++) {
2761
2762
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
2763
	for (int i = 0; i < nGrads; i++) {
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
2775
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}
	A = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void General_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2776
2777
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2778
2779
2780

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

Thomas Witkowski's avatar
Thomas Witkowski committed
2781
2782
    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
2783
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
2784
2785
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
2786
2787
2788
2789
    }
  }

  void General_FOT::getLb(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
2790
			  int nPoints, 
2791
2792
			  VectorOfFixVecs<DimVec<double> >& result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2793
2794
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2795

Thomas Witkowski's avatar
Thomas Witkowski committed
2796
2797
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2798
2799

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
2800

Thomas Witkowski's avatar
Thomas Witkowski committed
2801
2802
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
2803
2804
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2805
      for (int i = 0; i < nGrads; i++) {
2806
2807
2808
2809
2810
2811
2812
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lb(Lambda, (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg), 
	 result[iq], 1.0);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2813
2814
  void General_FOT::eval(int nPoints,
			 const double *uhAtQP,
2815
2816
2817
2818
2819
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
2820
    int dow = Global::getGeo(WORLD);  
Thomas Witkowski's avatar
Thomas Witkowski committed
2821
2822
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2823

Thomas Witkowski's avatar
Thomas Witkowski committed
2824
2825
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2826

2827
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2828
      for (int iq = 0; iq < nPoints; iq++) {
2829
2830
	double resultQP = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
2831
	for (int i = 0; i < nVecs; i++) {
2832
2833
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
2834
	for (int i = 0; i < nGrads; i++) {
2835
2836
2837
2838
2839
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}

	const WorldVector<double> &b = (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);

2840
	for (int i = 0; i < dow; i++) {
2841
2842
2843
2844
2845
2846
2847
2848
2849
2850
	  resultQP += grdUhAtQP[iq][i] * b[i];
	}
	result[iq] += factor * resultQP;
      }
    }
  }

  void General_ZOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
				Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2851
2852
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2853
2854
2855

    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

Thomas Witkowski's avatar
Thomas Witkowski committed
2856
    for (int i = 0; i < nVecs; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
2857
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
2858
    for (int i = 0; i < nGrads; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
2859
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
2860
2861
  }

2862
2863
  void General_ZOT::getC(const ElInfo *, int nPoints,
			 std::vector<double> &C) const
2864
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2865
2866
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2867

Thomas Witkowski's avatar
Thomas Witkowski committed
2868
2869
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2870

Thomas Witkowski's avatar
Thomas Witkowski committed
2871
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2872
      for (int i = 0; i < nVecs; i++)
2873
	vecsArg[i] = vecsAtQPs_[i][iq];
Thomas Witkowski's avatar
Thomas Witkowski committed
2874
2875
      for (int i = 0; i < nGrads; i++)

2876
	gradsArg[i] = gradsAtQPs_[i][iq];
Thomas Witkowski's avatar
Thomas Witkowski committed
2877

2878
2879
2880
2881
      C[iq] += (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2882
2883
  void General_ZOT::eval(int nPoints,
			 const double *uhAtQP,
2884
2885
2886
2887
2888
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2889
2890
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2891

Thomas Witkowski's avatar
Thomas Witkowski committed
2892
2893
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2894

Thomas Witkowski's avatar
Thomas Witkowski committed
2895
2896
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
2897
2898
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2899
      for (int i = 0; i < nGrads; i++) {
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      result[iq] += fac * (*f_)(coordsAtQPs_[iq], vecsArg, gradsArg) * uhAtQP[iq];
    }
  }



  //************************************ parametric general Operatorterms *****************
  void GeneralParametric_SOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2913
2914
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2915
2916
2917
2918

    elInfo->getElementNormal(elementNormal_);
    coordsAtQPs_ = subAssembler->getCoordsAtQPs(elInfo, quad);

Thomas Witkowski's avatar
Thomas Witkowski committed
2919
2920
    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
2921
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
2922
2923
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
2924
2925
2926
2927
    }
  }

  void GeneralParametric_SOT::getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
2928
				      int nPoints, 
2929
2930
				      DimMat<double> **LALt) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2931
2932
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2933

Thomas Witkowski's avatar
Thomas Witkowski committed
2934
2935
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2936
2937

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
2938

Thomas Witkowski's avatar
Thomas Witkowski committed
2939
2940
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
2941
2942
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2943
      for (int i = 0; i < nGrads; i++) {
2944
2945
2946
2947
2948
2949
2950
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      lalt(Lambda, (*f_)(coordsAtQPs_[iq], elementNormal_, vecsArg, gradsArg), 
	   *(LALt[iq]), symmetric_, 1.0);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
2951
2952
  void GeneralParametric_SOT::eval(int nPoints,
				   const double *uhAtQP,
2953
2954
2955
2956
2957
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double factor) const
  {
2958
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
2959
2960
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2961

Thomas Witkowski's avatar
Thomas Witkowski committed
2962
2963
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
2964

Thomas Witkowski's avatar
Thomas Witkowski committed
2965
    for (int iq = 0; iq < nPoints; iq++) {
2966
2967
      double resultQP = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
2968
      for (int i = 0; i < nVecs; i++) {
2969
2970
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2971
      for (int i = 0; i < nGrads; i++) {
2972
2973
2974
2975
2976
	gradsArg[i] = gradsAtQPs_[i][iq];
      }

      WorldMatrix<double> A = (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);

2977
2978
2979
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
2980
2981
2982
2983
2984
	    resultQP += A[i][j] * D2UhAtQP[iq][j][i];
	  }
	}
      }

2985
      if (grdUhAtQP) {
2986
2987
2988
2989
2990
2991
2992
	resultQP += (*divFct_)(A) * grdUhAtQP[iq];
      }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2993
  void GeneralParametric_SOT::weakEval(int nPoints,
2994
2995
2996
				       const WorldVector<double> *grdUhAtQP,
				       WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2997
2998
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2999

Thomas Witkowski's avatar
Thomas Witkowski committed
3000
    std::vector<double> vecsArg(nVecs);
For faster browsing, not all history is shown. View entire blame