Operator.cc 95.2 KB
Newer Older
3001

3002
    if (grdUhAtQP) {
3003
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
3004
3005
      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nVecs; i++) {
3006
3007
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
3008
	for (int i = 0; i < nGrads; i++) {
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}
	A = (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);
	result[iq] += A * grdUhAtQP[iq];
      }
    }
  }

  void GeneralParametric_FOT::initElement(const ElInfo* elInfo, 
					  SubAssembler* subAssembler,
					  Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3021
3022
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3023
3024
3025
3026

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3027
3028
    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
3029
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
3030
3031
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
3032
3033
3034
3035
    }
  }

  void GeneralParametric_FOT::getLb(const ElInfo *elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3036
				    int nPoints, 
3037
3038
				    VectorOfFixVecs<DimVec<double> >& result) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3039
3040
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3041

Thomas Witkowski's avatar
Thomas Witkowski committed
3042
3043
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
3044
3045

    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
3046
3047
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
3048
3049
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3050
      for (int i = 0; i < nGrads; i++) {
3051
3052
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
3053
      lb(Lambda, (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg), result[iq], 1.0);
3054
3055
3056
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3057
3058
  void GeneralParametric_FOT::eval(int nPoints,
				   const double *uhAtQP,
3059
3060
3061
3062
3063
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double factor) const
  {
3064
    int dow = Global::getGeo(WORLD);
Thomas Witkowski's avatar
Thomas Witkowski committed
3065
3066
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3067

Thomas Witkowski's avatar
Thomas Witkowski committed
3068
3069
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
3070

3071
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3072
      for (int iq = 0; iq < nPoints; iq++) {
3073
3074
	double resultQP = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
3075
	for (int i = 0; i < nVecs; i++) {
3076
3077
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
3078
	for (int i = 0; i < nGrads; i++) {
3079
3080
3081
3082
3083
	  gradsArg[i] = gradsAtQPs_[i][iq];
	}

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

3084
	for (int i = 0; i < dow; i++) {
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
	  resultQP += grdUhAtQP[iq][i] * b[i];
	}
	result[iq] += factor * resultQP;
      }
    }
  }

  void GeneralParametric_ZOT::initElement(const ElInfo* elInfo, SubAssembler* subAssembler,
					  Quadrature *quad)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3095
3096
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3097
3098
3099
3100

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3101
3102
    for (int i = 0; i < nVecs; i++) {
      vecsAtQPs_[i] = getVectorAtQPs(vecs_[i], elInfo, subAssembler, quad);
3103
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
3104
3105
    for (int i = 0; i < nGrads; i++) {
      gradsAtQPs_[i] = getGradientsAtQPs(grads_[i], elInfo, subAssembler, quad);
3106
3107
3108
    }
  }

3109
3110
  void GeneralParametric_ZOT::getC(const ElInfo *, int nPoints, 
				   std::vector<double> &C) const
3111
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3112
3113
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3114

Thomas Witkowski's avatar
Thomas Witkowski committed
3115
3116
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
3117

Thomas Witkowski's avatar
Thomas Witkowski committed
3118
3119
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
3120
3121
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3122
      for (int i = 0; i < nGrads; i++) {
3123
3124
3125
3126
3127
3128
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      C[iq] += (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg);
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3129
3130
  void GeneralParametric_ZOT::eval(int nPoints,
				   const double *uhAtQP,
3131
3132
3133
3134
3135
				   const WorldVector<double> *grdUhAtQP,
				   const WorldMatrix<double> *D2UhAtQP,
				   double *result,
				   double fac) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3136
3137
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3138

Thomas Witkowski's avatar
Thomas Witkowski committed
3139
3140
    std::vector<double> vecsArg(nVecs);
    std::vector<WorldVector<double> > gradsArg(nGrads);
3141

Thomas Witkowski's avatar
Thomas Witkowski committed
3142
3143
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
3144
3145
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3146
      for (int i = 0; i < nGrads; i++) {
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
	gradsArg[i] = gradsAtQPs_[i][iq];
      }
      result[iq] += fac * (*f_)(coordsAtQPs_[iq],  elementNormal_, vecsArg, gradsArg) * uhAtQP[iq];
    }
  }

  void VecAndVecOfGradientsAtQP_ZOT::initElement(const ElInfo* elInfo,
						 SubAssembler* subAssembler,
						 Quadrature *quad)
  {
3157
    int size = static_cast<int>(vecs.size());
3158
    for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
3159
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
3160

Thomas Witkowski's avatar
Thomas Witkowski committed
3161
    vecAtQPs = getVectorAtQPs(vec, elInfo, subAssembler, quad);
3162
3163
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3164
  void VecAndVecOfGradientsAtQP_ZOT::getC(const ElInfo *, int nPoints,
3165
					  std::vector<double> &C) const
3166
  {
3167
    int size = static_cast<int>(vecs.size());
3168
3169
    std::vector<WorldVector<double>*> arg(size);

Thomas Witkowski's avatar
Thomas Witkowski committed
3170
    for (int iq = 0; iq < nPoints; iq++) {
3171
      for (int i = 0; i < size; i++)
3172
	arg[i] = &(gradsAtQPs[i][iq]);
3173

3174
3175
      C[iq] += (*f)(vecAtQPs[iq], arg);
    }
3176
3177
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3178
3179
  void VecAndVecOfGradientsAtQP_ZOT::eval(int nPoints,
					  const double *uhAtQP,
3180
3181
3182
3183
3184
					  const WorldVector<double> *grdUhAtQP,
					  const WorldMatrix<double> *D2UhAtQP,
					  double *result,
					  double fac) const
  {
3185
    int size = static_cast<int>(vecs.size());
3186
3187
    std::vector<WorldVector<double>*> arg(size);

Thomas Witkowski's avatar
Thomas Witkowski committed
3188
    for (int iq = 0; iq < nPoints; iq++) {
3189
      for (int i = 0; i < size; i++)
3190
	arg[i] = &(gradsAtQPs[i][iq]);
3191

3192
3193
      result[iq] += fac * (*f)(vecAtQPs[iq], arg) * uhAtQP[iq];
    }
3194
3195
  }

3196
3197
3198
3199
  void VecGrad_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3200
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
    vecAtQPs = getVectorAtQPs(vec1, elInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec2, elInfo, subAssembler, quad);
  }

  void VecGrad_FOT::initElement(const ElInfo* smallElInfo,
				const ElInfo* largeElInfo,
				SubAssembler* subAssembler,
				Quadrature *quad)
  {
    vecAtQPs = getVectorAtQPs(vec1, smallElInfo, largeElInfo, subAssembler, quad);
    gradAtQPs = getGradientsAtQPs(vec2, smallElInfo, largeElInfo, subAssembler, quad);
3211
3212
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3213
3214
  void VecGrad_FOT::getLb(const ElInfo *elInfo, int nPoints, 
			  VectorOfFixVecs<DimVec<double> >& Lb) const {
3215
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
3216
    for (int iq = 0; iq < nPoints; iq++)
3217
3218
3219
      lb(Lambda, (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]), Lb[iq], 1.0);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
3220
  void VecGrad_FOT::eval(int nPoints,
3221
3222
3223
3224
3225
3226
3227
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3228
      for (int iq = 0; iq < nPoints; iq++) {
3229
3230
3231
3232
3233
3234
	WorldVector<double> b = (*vecFct)(vecAtQPs[iq], gradAtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3235
  void VecGrad_SOT::eval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
3236
3237
3238
3239
3240
3241
3242
3243
3244
				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
3245
      for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3246
3247
3248
3249
3250
3251
3252
3253
3254
3255
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);
	double resultQP = 0.0;
	for (int i = 0; i < dow; i++) {
	  resultQP += D2UhAtQP[iq][i][i];
	}
	result[iq] += fac * resultQP * factor;
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3256
  void VecGrad_SOT::weakEval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
3257
3258
3259
3260
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3261
      for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3262
3263
3264
3265
3266
3267
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3268
  void VecGrad_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
Thomas Witkowski's avatar
Thomas Witkowski committed
3269
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
3270
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3271
3272
3273
3274
3275
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
    }
  }
  
  void VecGrad_SOT::initElement(const ElInfo* elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3276
3277
				SubAssembler* subAssembler,
				Quadrature *quad) 
Thomas Witkowski's avatar
Thomas Witkowski committed
3278
3279
3280
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290
  {
    vecAtQPs = subAssembler->getVectorAtQPs(vec1, elInfo, quad);
    gradAtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
  }

  void FctGrad2_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
    grad1AtQPs = subAssembler->getGradientsAtQPs(vec1, elInfo, quad);
    grad2AtQPs = subAssembler->getGradientsAtQPs(vec2, elInfo, quad);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3291
  void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
Thomas Witkowski's avatar
Thomas Witkowski committed
3292
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
3293
    for (int iq = 0; iq < nPoints; iq++)
Thomas Witkowski's avatar
Thomas Witkowski committed
3294
3295
3296
      lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
3297
  void FctGrad2_FOT::eval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
3298
3299
3300
3301
3302
3303
3304
			 const double *uhAtQP,
			 const WorldVector<double> *grdUhAtQP,
			 const WorldMatrix<double> *D2UhAtQP,
			 double *result,
			 double factor) const
  {
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3305
      for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3306
3307
3308
3309
3310
3311
	WorldVector<double> b = (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

3312
}
For faster browsing, not all history is shown. View entire blame