DOFVector.hh 33.7 KB
Newer Older
1001
1002
1003
  {
    typename DOFVector<T>::Iterator vIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
1004
1005
1006
1007
1008
1009
    for(vIterator.reset(), rIterator.reset();
	!vIterator.end();
	++vIterator, ++rIterator) 
      {  
	*rIterator = (*vIterator) + scal;
      };
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
    return result;
  }

  template<typename T>
  inline const DOFVector<T>& add(const DOFVector<T>& v1,
				 const DOFVector<T>& v2,
				 DOFVector<T>& result)
  {
    typename DOFVector<T>::Iterator v1Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v1)), USED_DOFS);
    typename DOFVector<T>::Iterator v2Iterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&v2)), USED_DOFS);
    typename DOFVector<T>::Iterator rIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
1021
1022
1023
1024
1025
1026
    for(v1Iterator.reset(), v2Iterator.reset(), rIterator.reset();
	!v1Iterator.end();
	++v1Iterator, ++v2Iterator, ++rIterator) 
      {  
	*rIterator = (*v1Iterator) + (*v2Iterator);
      };
1027
    return result;
1028

1029
1030
1031
  }

  template<typename T>
1032
  const T *DOFVectorBase<T>::getLocalVector(const Element *el, T *d) const
1033
  {
1034
1035
    static T* localVec = NULL;
    const DOFAdmin* admin = feSpace->getAdmin();
1036

1037
1038
    T *result;
  
1039
    if (d) {
1040
1041
1042
1043
1044
1045
      result = d;
    } else {
      if(localVec) delete [] localVec;
      localVec = new T[nBasFcts]; 
      result = localVec;
    }
1046

1047
1048
    const DegreeOfFreedom *localIndices = 
      feSpace->getBasisFcts()->getLocalIndices(el, admin, NULL);
1049

1050
    for (int i = 0; i < nBasFcts; i++)
1051
      result[i] = (*this)[localIndices[i]];
1052
1053

    return result;
1054
1055
1056
  }

  template<typename T>
1057
1058
  const T *DOFVectorBase<T>::getVecAtQPs(const ElInfo         *elInfo, 
					 const Quadrature     *quad,
1059
					 const FastQuadrature *quadFast,
1060
					 T                    *vecAtQPs) const
1061
1062
1063
  {
    FUNCNAME("DOFVector<T>::getVecAtQPs()");
  
1064
    TEST_EXIT(quad || quadFast)("neither quad nor quadFast defined\n");
1065

1066
1067
1068
    if(quad && quadFast) {
      TEST_EXIT(quad == quadFast->getQuadrature())
	("quad != quadFast->quadrature\n");
1069
1070
    }

1071
    TEST_EXIT(!quadFast || quadFast->getBasisFunctions() == feSpace->getBasisFcts())
1072
1073
      ("invalid basis functions");

1074
1075
1076
1077
    Element *el = elInfo->getElement();

    const Quadrature *quadrature = quadFast ? quadFast->getQuadrature() : quad;

1078
    const BasisFunction *basFcts = feSpace->getBasisFcts();
1079
1080
1081
1082
1083

    int numPoints = quadrature->getNumPoints();
    int nBasFcts  = basFcts->getNumber();
    int i, j;

1084
    static T *localvec = NULL;
1085

1086
1087
1088
1089
1090
    T *result;

    if (vecAtQPs) {
      result = vecAtQPs;
    } else {
1091
1092
1093
      if(localvec) delete [] localvec;
      localvec = new T[numPoints];
      for(i = 0; i < numPoints; i++) {
1094
1095
1096
1097
	localvec[i] = 0.0;
      }
      result = localvec;
    }
1098
1099
1100
1101
1102
1103
  
    const T *localVec = getLocalVector(el, NULL);

    for(i = 0; i < numPoints; i++) {
      for(result[i] = j = 0; j < nBasFcts; j++) {
	result[i] +=  localVec[j] * 
1104
1105
1106
1107
1108
	  (quadFast ? 
	   (quadFast->getPhi(i, j)) : 
	   ((*(basFcts->getPhi(j)))(quad->getLambda(i))));
      }
    }
1109

1110
1111
1112
    return const_cast<const T*>(result);
  }

1113
1114
1115
1116
1117
1118
1119



  // 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
1120
  {
1121
1122
    return v.getSize();
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
1123

1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
  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
1143
1144
  }

1145
1146
1147
1148
1149
1150
  template<typename T>
  double DOFVector<T>::DoubleWell(Quadrature* q) const
  {
    FUNCNAME("DOFVector::DoubleWell()");
  
    if (!q) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1151
1152
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
1153
1154
1155
1156
1157
1158
    }

    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
1159
1160
1161
    Mesh* mesh = feSpace->getMesh();   
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, 
		   DoubleWell_fct);
1162
1163
1164
1165
1166
1167
1168

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::DoubleWell_fct(ElInfo *elinfo)
  {   
Thomas Witkowski's avatar
Thomas Witkowski committed
1169
1170
    const T *uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
    int nPoints = quad_fast->getNumPoints();
1171
    double normT = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
1172
1173
    for (int iq = 0; iq < nPoints; iq++) {
      normT += quad_fast->getWeight(iq) * sqr(uh_vec[iq]) * sqr(1.0 - uh_vec[iq]);
1174
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
1175
    norm += elinfo->getDet() * normT;
1176
1177
1178
1179
    
    return 0;
  }

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