Skip to content
Snippets Groups Projects
Commit 934bf0c1 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

pp

parent 83430f2b
No related branches found
No related tags found
No related merge requests found
...@@ -3,14 +3,15 @@ ...@@ -3,14 +3,15 @@
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>. // bugs and problems to <gmsh@geuz.org>.
#include "GmshMessage.h"
#include "GaussIntegration.h" #include "GaussIntegration.h"
#include "GaussLegendre1D.h" #include "GaussLegendre1D.h"
#include "stdio.h"
IntPt *getGQPriPts(int order); IntPt *getGQPriPts(int order);
int getNGQPriPts(int order); int getNGQPriPts(int order);
IntPt * GQP[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; IntPt * GQP[] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
IntPt *getGQPriPts(int order) IntPt *getGQPriPts(int order)
{ {
...@@ -18,13 +19,11 @@ IntPt *getGQPriPts(int order) ...@@ -18,13 +19,11 @@ IntPt *getGQPriPts(int order)
int nTri = getNGQTPts(order); int nTri = getNGQTPts(order);
int n = nLin*nTri; int n = nLin*nTri;
int index = order; int index = order;
if (order >= sizeof(GQP) / sizeof(IntPt*)) if (order >= (int)(sizeof(GQP) / sizeof(IntPt*)))
Msg::Fatal("Increase size of GQP in gauss quadrature prism"); Msg::Fatal("Increase size of GQP in gauss quadrature prism");
if(!GQP[index]) if(!GQP[index]){
{
double *linPt,*linWt; double *linPt,*linWt;
IntPt *triPts = getGQTPts(order); IntPt *triPts = getGQTPts(order);
// printf("o = %d n = %d nT = %d nL = %d i = %d\n",order,n,nTri,nLin,index);
gmshGaussLegendre1D(nLin,&linPt,&linWt); gmshGaussLegendre1D(nLin,&linPt,&linWt);
GQP[index] = new IntPt[n]; GQP[index] = new IntPt[n];
int l = 0; int l = 0;
...@@ -34,7 +33,6 @@ IntPt *getGQPriPts(int order) ...@@ -34,7 +33,6 @@ IntPt *getGQPriPts(int order)
GQP[index][l].pt[1] = triPts[i].pt[1]; GQP[index][l].pt[1] = triPts[i].pt[1];
GQP[index][l].pt[2] = linPt[j]; GQP[index][l].pt[2] = linPt[j];
GQP[index][l++].weight = triPts[i].weight*linWt[j]; GQP[index][l++].weight = triPts[i].weight*linWt[j];
// printf ("%d: %f %f %f %f\n",l-1,triPts[i].pt[0],triPts[i].pt[1],linPt[j],triPts[i].weight*linWt[j]);
} }
} }
} }
......
...@@ -542,7 +542,9 @@ class fullMatrix ...@@ -542,7 +542,9 @@ class fullMatrix
} }
// specific functions for dgshell // specific functions for dgshell
void mult_naiveBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol, const int alpha, const int beta, fullVector<scalar> &c, const int row=0) const void mult_naiveBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol,
const int alpha, const int beta, fullVector<scalar> &c,
const int row=0) const
{ {
if(beta != 1) if(beta != 1)
c.scale(beta); c.scale(beta);
...@@ -550,7 +552,8 @@ class fullMatrix ...@@ -550,7 +552,8 @@ class fullMatrix
for(int k = 0; k < _c ; k++) for(int k = 0; k < _c ; k++)
c._data[j] += alpha*(*this)(row, k) * b(k, j); c._data[j] += alpha*(*this)(row, k) * b(k, j);
} }
void multOnBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol, const int alpha, const int beta, fullVector<scalar> &c) const void multOnBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol,
const int alpha, const int beta, fullVector<scalar> &c) const
#if !defined(HAVE_BLAS) #if !defined(HAVE_BLAS)
{ {
mult_naiveBlock(b,ncol,fcol,alpha,beta,c); mult_naiveBlock(b,ncol,fcol,alpha,beta,c);
...@@ -558,7 +561,8 @@ class fullMatrix ...@@ -558,7 +561,8 @@ class fullMatrix
#endif #endif
; ;
void multWithATranspose(const fullVector<scalar> &x, scalar alpha, scalar beta, fullVector<scalar> &y) const void multWithATranspose(const fullVector<scalar> &x, scalar alpha, scalar beta,
fullVector<scalar> &y) const
#if !defined(HAVE_BLAS) #if !defined(HAVE_BLAS)
{ {
y.scale(beta); y.scale(beta);
...@@ -580,7 +584,8 @@ class fullMatrix ...@@ -580,7 +584,8 @@ class fullMatrix
scalar alpha=1., scalar beta=1.) scalar alpha=1., scalar beta=1.)
#if !defined(HAVE_BLAS) #if !defined(HAVE_BLAS)
{ {
Msg::Error("gemmWithAtranspose is only available with blas. If blas is not installed please transpose a before used gemm_naive"); Msg::Error("gemmWithAtranspose is only available with blas. If blas is not "
"installed please transpose a before used gemm_naive");
} }
#endif #endif
; ;
......
...@@ -82,7 +82,8 @@ the <a href="/gmsh/doc/LICENSE.txt">GNU General Public License ...@@ -82,7 +82,8 @@ the <a href="/gmsh/doc/LICENSE.txt">GNU General Public License
<li>Development version: <li>Development version:
<ul><li>automated nightly builds <ul><li>automated nightly builds
(<a href="http://onelab.info/CDash/index.php?project=Gmsh">dashboard</a>): (<a href="http://onelab.info/CDash/index.php?project=Gmsh">dashboard</a>):
<a href="/gmsh/bin/Windows/gmsh-svn-Windows.zip">Windows</a>, Windows <a href="/gmsh/bin/Windows/gmsh-svn-Windows.zip">32 bit</a> /
<a href="/gmsh/bin/Windows/gmsh-svn-Windows64.zip">64 bit</a>,
Linux <a href="/gmsh/bin/Linux/gmsh-svn-Linux.tgz">32 bit</a> / Linux <a href="/gmsh/bin/Linux/gmsh-svn-Linux.tgz">32 bit</a> /
<a href="/gmsh/bin/Linux/gmsh-svn-Linux64.tgz">64 bit</a> <a href="/gmsh/bin/Linux/gmsh-svn-Linux64.tgz">64 bit</a>
and and
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment