Skip to content

Commit

Permalink
Removed wrong scaling of Gauss weights
Browse files Browse the repository at this point in the history
  • Loading branch information
iBatistic committed Jun 18, 2024
1 parent 419d11e commit e81804c
Showing 1 changed file with 4 additions and 7 deletions.
11 changes: 4 additions & 7 deletions src/foam/primitives.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@
Description
Different geometrical primitives
"""
import sys

import numpy as np

Expand Down Expand Up @@ -60,16 +61,16 @@ def GaussPointsAndWeights(self, GaussPoinsNb, emptyDir) -> ([np.ndarray], [np.nd
if emptyDir is None:
# 3D case, Gauss points are distributed on face
# TO_DO
pass
sys.exit(1)
else:
# 2D case, Gauss points are distributed on the line
GaussPoints, weights = np.polynomial.legendre.leggauss(GaussPoinsNb)

facePoints = self.points()

# Normalize the weights
weights /= np.sum(weights)

facePoints = self.points()

# Project faces onto normal plane defined with emtyDir
faceIn2D = [point[:emptyDir] for point in facePoints]

Expand All @@ -89,9 +90,7 @@ def GaussPointsAndWeights(self, GaussPoinsNb, emptyDir) -> ([np.ndarray], [np.nd
faceIn2D = [list(inner_tuple) for inner_tuple in unique_lists]

# Get position of gauss points using face points
# Scale weights of the quadrature because we are not integrating from -1 to 1
faceGaussPoints = []
faceWeights = []

for gp,w in zip(GaussPoints, weights):
pointA = np.array(faceIn2D[1])
Expand All @@ -100,10 +99,8 @@ def GaussPointsAndWeights(self, GaussPoinsNb, emptyDir) -> ([np.ndarray], [np.nd
halfFacePoint = pointB + (pointA - pointB)/2

faceGaussPoint = halfFacePoint + halfFaceLen*gp*(pointA-pointB)/(2*halfFaceLen)
faceWeight = w * halfFaceLen

faceGaussPoints.append(faceGaussPoint)
faceWeights.append(faceWeight)

# Face centre in empty direction
faceCentreEmptyDir = self.centre()[emptyDir]
Expand Down

0 comments on commit e81804c

Please sign in to comment.