Operator.cc 95.2 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
3001
    std::vector<WorldVector<double> > gradsArg(nGrads);
3002

3003
    if (grdUhAtQP) {
3004
      WorldMatrix<double> A;
Thomas Witkowski's avatar
Thomas Witkowski committed
3005
3006
      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nVecs; i++) {
3007
3008
	  vecsArg[i] = vecsAtQPs_[i][iq];
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
3009
	for (int i = 0; i < nGrads; i++) {
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
	  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
3022
3023
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3024
3025
3026
3027

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

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

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

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

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

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

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

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

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

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

3085
	for (int i = 0; i < dow; i++) {
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
	  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
3096
3097
    int nVecs = static_cast<int>(vecs_.size());
    int nGrads = static_cast<int>(grads_.size());
3098
3099
3100
3101

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3143
3144
    for (int iq = 0; iq < nPoints; iq++) {
      for (int i = 0; i < nVecs; i++) {
3145
3146
	vecsArg[i] = vecsAtQPs_[i][iq];
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
3147
      for (int i = 0; i < nGrads; i++) {
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
	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)
  {
3158
    int size = static_cast<int>(vecs.size());
3159
    for (int i = 0; i < size; i++)
Thomas Witkowski's avatar
Thomas Witkowski committed
3160
      gradsAtQPs[i] = getGradientsAtQPs(vecs[i], elInfo, subAssembler, quad);
3161

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

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

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

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

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

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

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

3197
3198
3199
3200
  void VecGrad_FOT::initElement(const ElInfo* elInfo, 
				SubAssembler* subAssembler,
				Quadrature *quad) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
3201
3202
3203
3204
3205
3206
3207
3208
3209
3210
3211
    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);
3212
3213
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
3236
  void VecGrad_SOT::eval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
3237
3238
3239
3240
3241
3242
3243
3244
3245
				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
3246
      for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3247
3248
3249
3250
3251
3252
3253
3254
3255
3256
	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
3257
  void VecGrad_SOT::weakEval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
3258
3259
3260
3261
				    const WorldVector<double> *grdUhAtQP,
				    WorldVector<double> *result) const
  {
    if (grdUhAtQP) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3262
      for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3263
3264
3265
3266
3267
3268
	double factor = (*f)(vecAtQPs[iq], gradAtQPs[iq]);  
	axpy(factor, grdUhAtQP[iq], result[iq]);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
3269
  void VecGrad_SOT::getLALt(const ElInfo *elInfo, int nPoints, DimMat<double> **LALt) const {
Thomas Witkowski's avatar
Thomas Witkowski committed
3270
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
3271
    for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3272
3273
3274
3275
3276
      l1lt(Lambda, *(LALt[iq]), (*f)(vecAtQPs[iq], gradAtQPs[iq]));
    }
  }
  
  void VecGrad_SOT::initElement(const ElInfo* elInfo, 
Thomas Witkowski's avatar
Thomas Witkowski committed
3277
3278
				SubAssembler* subAssembler,
				Quadrature *quad) 
Thomas Witkowski's avatar
Thomas Witkowski committed
3279
3280
3281
3282
3283
3284
3285
3286
3287
3288
3289
3290
3291
  {
    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
3292
  void FctGrad2_FOT::getLb(const ElInfo *elInfo, int nPoints, VectorOfFixVecs<DimVec<double> >& Lb) const {
Thomas Witkowski's avatar
Thomas Witkowski committed
3293
    const DimVec<WorldVector<double> > &Lambda = elInfo->getGrdLambda();
Thomas Witkowski's avatar
Thomas Witkowski committed
3294
    for (int iq = 0; iq < nPoints; iq++)
Thomas Witkowski's avatar
Thomas Witkowski committed
3295
3296
3297
      lb(Lambda, (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]), Lb[iq], 1.0);
  }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
3298
  void FctGrad2_FOT::eval(int nPoints,
Thomas Witkowski's avatar
Thomas Witkowski committed
3299
3300
3301
3302
3303
3304
3305
			 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
3306
      for (int iq = 0; iq < nPoints; iq++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
3307
3308
3309
3310
3311
3312
	WorldVector<double> b = (*vecFct)(grad1AtQPs[iq], grad2AtQPs[iq]);
	result[iq] += b * grdUhAtQP[iq] * factor;
      }
    }
  }

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