DOFVector.hh 31.3 KB
Newer Older
1001
1002
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
1003
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
1004

Thomas Witkowski's avatar
Thomas Witkowski committed
1005
    if (quad && quadFast)
1006
1007
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
1008

1009
    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1010
1011
      ("invalid basis functions");

1012
1013
    Element *el = elInfo->getElement();
    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;
1014
    const BasisFunction *basFcts = feSpace->getBasisFcts();
1015
1016
    int numPoints = quadrature->getNumPoints();
    int nBasFcts  = basFcts->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
1017
    int j;
1018
1019
1020
1021
1022
1023
    static T *localvec = NULL;
    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
1024
1025
      if (localvec) 
	delete [] localvec;
1026
      localvec = new T[numPoints];
Thomas Witkowski's avatar
Thomas Witkowski committed
1027
      for (int i = 0; i < numPoints; i++)
1028
	localvec[i] = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1029

1030
1031
      result = localvec;
    }
1032
1033
1034
  
    const T *localVec = getLocalVector(el, NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
1035
1036
1037
    for (int i = 0; i < numPoints; i++) {
      for (result[i] = j = 0; j < nBasFcts; j++) {
	result[i] += localVec[j] * 
1038
1039
1040
1041
1042
	  (quadFast ? 
	   (quadFast->getPhi(i, j)) : 
	   ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
      }
    }
1043

1044
1045
1046
    return const_cast<const T*>(result);
  }

1047
1048
1049
1050
1051
1052
1053



  // Some free functions used in MTL4

  template <typename T>
  inline std::size_t size(const AMDiS::DOFVector<T>& v)
Thomas Witkowski's avatar
Thomas Witkowski committed
1054
  {
1055
1056
    return v.getSize();
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1057

1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
  template <typename T>
  inline std::size_t num_rows(const AMDiS::DOFVector<T>& v)
  {
    return v.getSize();
  }

  template <typename T>
  inline std::size_t num_cols(const AMDiS::DOFVector<T>& v)
  {
    return 1;
  }
 
  template <typename T>
  inline void set_to_zero(AMDiS::DOFVector<T>& v)
  {
    using math::zero;
    T  ref, my_zero(zero(ref));

    std::fill(v.begin(), v.end(), my_zero);
Thomas Witkowski's avatar
Thomas Witkowski committed
1077
1078
  }

1079
1080
1081
1082
1083
1084
  template<typename T>
  double DOFVector<T>::DoubleWell(Quadrature* q) const
  {
    FUNCNAME("DOFVector::DoubleWell()");
  
    if (!q) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1085
1086
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
1087
1088
1089
1090
1091
1092
    }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
1093
1094
1095
    Mesh* mesh = feSpace->getMesh();   
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   DoubleWell_fct);
1096
1097
1098
1099
1100
1101
1102

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
  {   
Thomas Witkowski's avatar
Thomas Witkowski committed
1103
1104
    const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int nPoints = quad_fast->getNumPoints();
1105
    double normT = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1106
    for (int iq = 0; iq < nPoints; iq++)
Thomas Witkowski's avatar
Thomas Witkowski committed
1107
      normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);
Thomas Witkowski's avatar
Thomas Witkowski committed
1108

Thomas Witkowski's avatar
Thomas Witkowski committed
1109
    norm += elinfo->getDet() * normT;
1110
1111
1112
1113
    
    return 0;
  }

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