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

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2114
  void VecMatrixGradientAtQP_SOT::eval(int nPoints,
2115
2116
2117
2118
2119
2120
2121
2122
				       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
2123
    for (int iq = 0; iq < nPoints; iq++) {
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
      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;
    }
  }
2143

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2204
2205
2206
      if (D2UhAtQP) {
	for (int i = 0; i < dow; i++) {
	  for (int j = 0; j < dow; j++) {
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
	    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
2218
  }
2219

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

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

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

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

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

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

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

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

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

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

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

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

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

  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
2408
				  int nPoints, 
2409
2410
2411
				  DimMat<double> **LALt) const
  {
    const DimVec<WorldVector<double> > Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
2412
    for (int iq = 0; iq < nPoints; iq++) {
2413
2414
2415
2416
      lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*g)(coordsAtQPs[iq]));
    }
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2432
  void CoordsAtQP_IJ_SOT::weakEval(int nPoints,
2433
2434
2435
				   const WorldVector<double> *grdUhAtQP,
				   WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2436
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2437
      for (int iq = 0; iq < nPoints; iq++) {
2438
2439
2440
2441
2442
2443
2444
2445
2446
2447
	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
2448
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
2449
2450
2451
  }
 
  void VecAtQP_IJ_SOT::getLALt(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
2452
			       int nPoints, 
2453
2454
2455
			       DimMat<double> **LALt) const
  {
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
2456
    for (int iq = 0; iq < nPoints; iq++) {
2457
2458
2459
2460
      lalt_kl(Lambda, xi, xj, *(LALt[iq]), (*f)(vecAtQPs[iq]));
    }
  }  

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2476
  void VecAtQP_IJ_SOT::weakEval(int nPoints,
2477
2478
2479
				const WorldVector<double> *grdUhAtQP,
				WorldVector<double> *result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2480
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2481
      for (int iq = 0; iq < nPoints; iq++) {
2482
2483
2484
2485
2486
2487
2488
2489
2490
2491
2492
	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
2493
2494
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2495
      vecsAtQPs[i] = getVectorAtQPs(vecs[i], elInfo, subAssembler, quad);
2496
2497
2498
    }
  }
 
2499
2500
  void VecOfDOFVecsAtQP_ZOT::getC(const ElInfo *, int nPoints,
				  std::vector<double> &C) const 
2501
  { 
2502
    int size = static_cast<int>(vecs.size());
2503
    std::vector<double> arg(size);
2504

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2523
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2524
      for (int i = 0; i < size; i++) {
2525
2526
2527
2528
2529
2530
2531
2532
2533
2534
2535
2536
2537
	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
2538
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
2539
2540
2541
2542
2543
2544
2545
    }
  }
 
  void VecDivergence_ZOT::initElement(const ElInfo* elInfo, 
				      SubAssembler* subAssembler,
				      Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
2546
2547
    int size = static_cast<int>(vecs.size());
    for (int i = 0; i < size; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
2548
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
2549
2550
2551
2552
    }
  }
 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2756
    if (grdUhAtQP) {
2757
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
2758
2759
      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nVecs; i++) {
2760
2761
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
2762
	for (int i = 0; i < nGrads; i++) {
2763
2764
2765
2766
2767
2768
2769
2770
2771
2772
2773
2774
	  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
2775
2776
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2777
2778
2779

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

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

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

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

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

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

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

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

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

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

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

2839
	for (int i = 0; i < dow; i++) {
2840
2841
2842
2843
2844
2845
2846
2847
2848
2849
	  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
2850
2851
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2852
2853
2854

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

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

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2894
2895
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
2896
2897
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2898
      for (int i = 0; i < nGrads; i++) {
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
	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
2912
2913
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2914
2915
2916
2917

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
2938
2939
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
2940
2941
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
2942
      for (int i = 0; i < nGrads; i++) {
2943
2944
2945
2946
2947
2948
2949
	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
2950
2951
  void GeneralParametric_SOT::eval(int nPoints,
				   const double *uhAtQP,
2952
2953
2954
2955
2956
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double factor) const
  {
2957
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
2958
2959
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
2960

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

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

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

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

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

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

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

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

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