/[ff3d]/ff3d/language/FunctionExpression.cpp
ViewVC logotype

Diff of /ff3d/language/FunctionExpression.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by delpinux, Thu Feb 20 19:32:09 2003 UTC revision 1.3 by delpinux, Fri Feb 21 13:34:07 2003 UTC
# Line 163  const real_t FunctionExpressionFEM::valu Line 163  const real_t FunctionExpressionFEM::valu
163    }    }
164  }  }
165    
166    
167    
168  template <typename MeshType>  template <typename MeshType>
169  const real_t FunctionExpressionFEM::  const real_t FunctionExpressionFEM::
170  __dx(const double& x,  __dx(const double& x,
# Line 180  __dx(const double& x, Line 182  __dx(const double& x,
182    
183    typedef typename MeshType::ElementsGeometry CellType;    typedef typename MeshType::ElementsGeometry CellType;
184    typedef typename FiniteElementTraits<CellType>::Transformation TransformationType;    typedef typename FiniteElementTraits<CellType>::Transformation TransformationType;
   typedef typename FiniteElementTraits<CellType>::JacobianTransformation Jacobian;  
185    typedef typename FiniteElementTraits<CellType>::Type FEMType;    typedef typename FiniteElementTraits<CellType>::Type FEMType;
186    
187    TransformationType T(*i);    TransformationType T(*i);
188    Jacobian J(T);  
189      TinyVector<3,real_t> Xhat;
190      if (not(T.invertT(x,y,z,Xhat))) {
191        return 0;
192      }
193    
194      TinyVector<CellType::NumberOfVertices, real_t> fiValues;
195    
196      TinyMatrix<3, CellType::NumberOfVertices> DwValues;
197    
198      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {
199        DwValues(0,k) = FEMType::instance().dxW(k,Xhat);
200        DwValues(1,k) = FEMType::instance().dyW(k,Xhat);
201        DwValues(2,k) = FEMType::instance().dzW(k,Xhat);
202      }
203    
204      Vector<real_t>& f = (*__value);
205      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {
206        fiValues[k] = f[M.vertexNumber((*i)(k))];
207      }
208    
209      TinyVector<3, real_t> gradientPart = DwValues*fiValues;
210    
211      TinyMatrix<3,3, real_t> J;
212      TinyVector<3, real_t> temp;
213    
214      T.dx(x,y,z,temp);
215      for(size_t i=0; i<3; ++i) {
216        J(0,i) = temp[i];
217      }
218    
219      T.dy(x,y,z,temp);
220      for(size_t i=0; i<3; ++i) {
221        J(1,i) = temp[i];
222      }
223    
224      T.dz(x,y,z,temp);
225      for(size_t i=0; i<3; ++i) {
226        J(2,i) = temp[i];
227      }
228    
229      // now we use
230    
231      TinyVector<3, real_t> result;
232    
233      gaussPivot(J, gradientPart, result);
234    
235      switch (d) {
236      case FunctionExpressionDx::x: {
237        return result[0];
238      }
239      case FunctionExpressionDx::y: {
240        return result[1];
241      }
242      case FunctionExpressionDx::z: {
243        return result[2];
244      }
245      default: {
246        fferr(0) << __FILE__ << ':' << __LINE__ << ": "
247                 << "trying to compute unknown derivative\n";
248        std::exit(1);
249        return 0;
250      }
251      }
252    }
253    
254    template <>
255    const real_t FunctionExpressionFEM::
256    __dx<Structured3DMesh>(const double& x,
257                           const double& y,
258                           const double& z,
259                           const FunctionExpressionDx::Direction& d) const
260    {
261      Structured3DMesh& M = static_cast<Structured3DMesh&>(*(*__mesh).mesh());
262    
263      Structured3DMesh::iterator i = M.find(x,y,z);
264    
265      if (i == 0) { // outside of the mesh!
266        return 0;
267      }
268    
269      typedef Structured3DMesh::ElementsGeometry CellType;
270      typedef FiniteElementTraits<CellType>::Transformation TransformationType;
271      typedef FiniteElementTraits<CellType>::Type FEMType;
272    
273      TransformationType T(*i);
274    
275    TinyVector<3,real_t> Xhat;    TinyVector<3,real_t> Xhat;
276    if (not(T.invertT(x,y,z,Xhat))) {    if (not(T.invertT(x,y,z,Xhat))) {
# Line 199  __dx(const double& x, Line 285  __dx(const double& x,
285      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {
286        DwiValues[k] = FEMType::instance().dxW(k,Xhat);        DwiValues[k] = FEMType::instance().dxW(k,Xhat);
287      }      }
288        TinyVector<3, real_t> Dx;
289        T.dx(x,y,z, Dx);
290        DwiValues /= Dx[0];
291      break;      break;
292    }    }
293    case FunctionExpressionDx::y: {    case FunctionExpressionDx::y: {
294      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {
295        DwiValues[k] = FEMType::instance().dyW(k,Xhat);        DwiValues[k] = FEMType::instance().dyW(k,Xhat);
296      }      }
297        TinyVector<3, real_t> Dy;
298        T.dy(x,y,z, Dy);
299        DwiValues /= Dy[1];
300      break;      break;
301    }    }
302    case FunctionExpressionDx::z: {    case FunctionExpressionDx::z: {
303      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {      for (size_t k = 0; k<CellType::NumberOfVertices; ++k) {
304        DwiValues[k] = FEMType::instance().dzW(k,Xhat);        DwiValues[k] = FEMType::instance().dzW(k,Xhat);
305      }      }
306        TinyVector<3, real_t> Dz;
307        T.dz(x,y,z, Dz);
308        DwiValues /= Dz[2];
309      break;      break;
310    }    }
311    default: {    default: {
# Line 225  __dx(const double& x, Line 320  __dx(const double& x,
320      fiValues[k] = f[M.vertexNumber((*i)(k))];      fiValues[k] = f[M.vertexNumber((*i)(k))];
321    }    }
322    
323    return (DwiValues*fiValues/J.jacobian(0,0));    return (DwiValues*fiValues);
324  }  }
325    
326  const real_t FunctionExpressionFEM::dx(const double& x,  const real_t FunctionExpressionFEM::dx(const double& x,
# Line 238  const real_t FunctionExpressionFEM::dx(c Line 333  const real_t FunctionExpressionFEM::dx(c
333      return __dx<Structured3DMesh>(x,y,z,d);      return __dx<Structured3DMesh>(x,y,z,d);
334      break;      break;
335    }    }
336      case Mesh::hexahedraMesh: {
337        return __dx<MeshOfHexahedra>(x,y,z,d);
338        break;
339      }
340    default: {    default: {
341      fferr(0) << __FILE__ << ':' << __LINE__      fferr(0) << __FILE__ << ':' << __LINE__
342               << ": not implemented\n";               << ": not implemented\n";

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26