diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp index f1444dc57d967d2da5041219062bd89e7cc77a3d..f2a0c8c85f8b29a11819eaf1f8b967eef95e2b29 100644 --- a/Geo/GFaceCompound.cpp +++ b/Geo/GFaceCompound.cpp @@ -901,7 +901,7 @@ void GFaceCompound::compute_distance() const double mu = L / 28; simpleFunction<double> DIFF(mu * mu), ONE(1.0); dofManager<double, double> myAssembler(_lsys); - distanceTerm distance(model(), 1, 1, &DIFF, &ONE); + distanceTerm distance(model(), 1, &DIFF, &ONE); std::vector<MVertex*> ordered; boundVertices(_U0, ordered); diff --git a/Solver/distanceTerm.h b/Solver/distanceTerm.h new file mode 100644 index 0000000000000000000000000000000000000000..4843889a515e0b413e97853d3c5228410943b024 --- /dev/null +++ b/Solver/distanceTerm.h @@ -0,0 +1,43 @@ +// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle +// +// See the LICENSE.txt file for license information. Please report all +// bugs and problems to <gmsh@geuz.org>. + +#ifndef _DISTANCE_TERM_H_ +#define _DISTANCE_TERM_H_ + +#include "helmholtzTerm.h" + +class distanceTerm : public helmholtzTerm<double> { + public: + distanceTerm(GModel *gm, int iFieldR, int iFieldC, + simpleFunction<double> *k, simpleFunction<double> *a) + : helmholtzTerm<double>(gm, iFieldR, iFieldC, k, a) {} + void elementVector(SElement *se, fullVector<double> &m) const + { + MElement *e = se->getMeshElement(); + int nbNodes = e->getNumVertices(); + int integrationOrder = 2 * e->getPolynomialOrder(); + int npts; + IntPt *GP; + double jac[3][3]; + double ff[256]; + e->getIntegrationPoints(integrationOrder, &npts, &GP); + + m.scale(0.); + + for (int i = 0; i < npts; i++){ + const double u = GP[i].pt[0]; + const double v = GP[i].pt[1]; + const double w = GP[i].pt[2]; + const double weight = GP[i].weight; + const double detJ = e->getJacobian(u, v, w, jac); + e->getShapeFunctions(u, v, w, ff); + for (int j = 0; j < nbNodes; j++){ + m(j) += ff[j] * weight * detJ; + } + } + } +}; + +#endif