From 453efcfad3d70b7bb653efbae1cb9a37eba3ba11 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Fri, 18 Jun 2010 10:51:14 +0000
Subject: [PATCH]
---
CMakeLists.txt | 2 +-
Solver/crossConfTerm.h | 4 +--
Solver/helmholtzTerm.h | 2 +-
Solver/orthogonalTerm.h | 75 +++++++++++++++++++++++++++++++++++++++++
4 files changed, 79 insertions(+), 4 deletions(-)
create mode 100644 Solver/orthogonalTerm.h
diff --git a/CMakeLists.txt b/CMakeLists.txt
index ac7414c412..ab89bf5dc8 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -79,7 +79,7 @@ set(GMSH_API
Geo/Homology.h
Mesh/meshGEdge.h Mesh/meshGFace.h Mesh/meshGFaceOptimize.h
Mesh/meshGFaceDelaunayInsertion.h
- Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h
+ Solver/dofManager.h Solver/femTerm.h Solver/laplaceTerm.h Solver/elasticityTerm.h Solver/crossConfTerm.h Solver/orthogonalTerm.h
Solver/linearSystem.h Solver/linearSystemGMM.h Solver/linearSystemCSR.h
Solver/linearSystemFull.h Solver/elasticitySolver.h
Post/PView.h Post/PViewData.h Plugin/PluginManager.h
diff --git a/Solver/crossConfTerm.h b/Solver/crossConfTerm.h
index f0dc4ddd80..a4a2146bde 100644
--- a/Solver/crossConfTerm.h
+++ b/Solver/crossConfTerm.h
@@ -46,7 +46,7 @@ class crossConfTerm : public femTerm<double> {
{
MElement *e = se->getMeshElement();
int nbNodes = e->getNumVertices();
- int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
+ int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
int npts;
IntPt *GP;
double jac[3][3];
@@ -54,7 +54,7 @@ class crossConfTerm : public femTerm<double> {
SVector3 Grads [256];
double grads[256][3];
e->getIntegrationPoints(integrationOrder, &npts, &GP);
-
+
m.setAll(0.);
for (int i = 0; i < npts; i++){
diff --git a/Solver/helmholtzTerm.h b/Solver/helmholtzTerm.h
index 2af2587d62..64b6d1890b 100644
--- a/Solver/helmholtzTerm.h
+++ b/Solver/helmholtzTerm.h
@@ -15,7 +15,7 @@
#include "fullMatrix.h"
#include "Numeric.h"
-// \nabla \cdot k \nabla U - a U = 0
+// \nabla \cdot k \nabla U - a U
template<class scalar>
class helmholtzTerm : public femTerm<scalar> {
protected:
diff --git a/Solver/orthogonalTerm.h b/Solver/orthogonalTerm.h
new file mode 100644
index 0000000000..cf1d3c3274
--- /dev/null
+++ b/Solver/orthogonalTerm.h
@@ -0,0 +1,75 @@
+// Gmsh - Copyright (C) 1997-2010 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _ORTHOGONAL_TERM_H_
+#define _ORTHOGONAL_TERM_H_
+
+#include "helmholtzTerm.h"
+
+class orthogonalTerm : public helmholtzTerm<double> {
+ protected:
+ fullVector<double> *_value;
+ public:
+ orthogonalTerm(GModel *gm, int iField, fullVector<double> &value)
+ : helmholtzTerm<double>(gm, iField, iField, 1.0, 0), _value(value) {}
+ void elementVector(SElement *se, fullVector<double> &m) const
+ {
+
+ MElement *e = se->getMeshElement();
+
+ //fill elementary matrix mat(i,j)
+ int nbNodes = e->getNumVertices();
+ int integrationOrder = 2 * (e->getPolynomialOrder() - 1);
+ int npts;
+ IntPt *GP;
+ double jac[3][3];
+ double invjac[3][3];
+ SVector3 Grads [256];
+ double grads[256][3];
+ e->getIntegrationPoints(integrationOrder, &npts, &GP);
+ fullMatrix<double> mat;
+ mat.setAll(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);
+ SPoint3 p; e->pnt(u, v, w, p);
+ const double _diff = (*_diffusivity)(p.x(), p.y(), p.z());
+ inv3x3(jac, invjac);
+ e->getGradShapeFunctions(u, v, w, grads);
+ for (int j = 0; j < nbNodes; j++){
+ Grads[j] = SVector3(invjac[0][0] * grads[j][0] + invjac[0][1] * grads[j][1] +
+ invjac[0][2] * grads[j][2],
+ invjac[1][0] * grads[j][0] + invjac[1][1] * grads[j][1] +
+ invjac[1][2] * grads[j][2],
+ invjac[2][0] * grads[j][0] + invjac[2][1] * grads[j][1] +
+ invjac[2][2] * grads[j][2]);
+ }
+ SVector3 N (jac[2][0], jac[2][1], jac[2][2]);
+ for (int j = 0; j < nbNodes; j++){
+ for (int k = 0; k <= j; k++){
+ mat(j, k) += dot(crossprod(Grads[j], Grads[k]), N) * weight * detJ * _diff;
+ }
+ }
+ }
+ for (int j = 0; j < nbNodes; j++)
+ for (int k = 0; k < j; k++)
+ mat(k, j) = -1.* m(j, k);
+
+ //2) compute vector m(i) = mat(i,j)*val(j)
+ fullVector<double> val(nbNodes);
+
+ m.scale(0.);
+ for (int i = 0; i < nbNodes; i++)
+ for (int j = 0; j < nbNodes; j++)
+ m(i) += mat(i,j)*val(j);
+
+ }
+};
+
+#endif
--
GitLab