<commit_before>from sys import platform
from os import path
from ctypes import CDLL, c_double, c_void_p, c_int, c_bool
import mozart as mz
import numpy as np

# OS Detection Code
prefix = "linux"
if platform == "linux" or platform == "linux32":
	prefix = "linux"
elif platform == "darwin":
	prefix = "osx"
elif platform == "win32":
	prefix = "win64"

dllpath = path.join(mz.__path__[0], prefix + '_' + 'libmozart.so')
lib = CDLL(dllpath)

def getMatrix1D(degree):
	"""
	Get FEM matrices on the reference domain

	Paramters
		- ``degree`` : degree of polynomial

	Returns
		- ``M_R`` : 
		- ``S_R`` : 
		- ``D_R`` :  

	"""
	M_R = None
	S_R = None
	D_R = None
	if degree is 1:
		M_R = np.array([[ 2,  1], [ 1, 2]]) / 3.
		S_R = np.array([[ 1, -1], [-1, 1]]) / 2.
		D_R = np.array([[-1,  1], [-1, 1]]) / 2.
	if degree is 2:
		M_R = np.array([[ 4,  2, -1], [ 2, 16,  2], [-1,  2, 4]]) / 15.
		S_R = np.array([[ 7, -8,  1], [-8, 16, -8], [ 1, -8, 7]]) / 6.
		D_R = np.array([[-3,  4, -1], [-1,  0,  1], [ 1, -4, 3]]) / 2.

	return (M_R, S_R, D_R)

def nJacobiP(x=np.array([0]), alpha=0, beta=0, degree=0):
	"""
	Evaluate normalized Jacobi polynomial of type alpha, beta > -1 at point x for order n
	(the special case of alpha = beta = 0, knows as the normalized Legendre polynomial)

	Paramters
		- ``x`` (``float64 array``) : variable x
		- ``alpha`` (``int32``) : superscript alpha of normalized Jacobi polynomial
		- ``beta`` (``int32``) : superscript beta of normalized Jacobi polynomial
		- ``degree`` (``int32``) : Polynomial degree

	Returns
		- ``P`` (``float64 array``) : the value of degree-th order normalized Jacobi polynomial at x
	
	Example
		>>> N = 2
		>>> x = np.array([-1, 0, 1])
		>>> from mozart.poisson.solve import nJacobiP
		>>> P = nJacobiP(x,0,0,N)
		>>> print(P)
		array([ 1.58113883, -0.79056942,  1.58113883])	
	"""
	Pn = np.zeros((degree+1,x.size),float)
	Pn[0,:] = np.sqrt(2**(-alpha-beta-1) * np.math.gamma(alpha+beta+2) / ((np.math.gamma(alpha + 1) * np.math.gamma(beta + 1))))

	if degree == 0:
		P = Pn
	else:
		Pn[1,:] = np.multiply(Pn[0,:]*np.sqrt((alpha+beta+3.0)/((alpha+1)*(beta+1))),((alpha+beta+2)*x+(alpha-beta)))/2
		a_n = 2.0/(2+alpha+beta)*np.sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
		for n in range(2,degree+1):
			anew=2.0/(2*n+alpha+beta)*np.sqrt(n*(n+alpha+beta)*(n+alpha)*(n+beta)/((2*n+alpha+beta-1)*(2*n+alpha+beta+1)))
			b_n=-(alpha**2-beta**2)/((2*(n-1)+alpha+beta)*(2*(n-1)+alpha+beta+2.0))
			Pn[n,:]=(np.multiply((x-b_n),Pn[n-1,:])-a_n*Pn[n-2,:])/anew;
			a_n=anew;

	P = Pn[degree,:]
	return P

def one_dim(c4n, n4e, n4Db, f, u_D, degree = 1):
	"""
	Computes the coordinates of nodes and elements.
	
	Parameters
		- ``c4n`` (``float64 array``) : coordinates
		- ``n4e`` (``int32 array``) : nodes for elements<commit_msg>function DnJacobiP is added<commit_after>from sys import platform
from os import path
from ctypes import CDLL, c_double, c_void_p, c_int, c_bool
import mozart as mz
import numpy as np

# OS Detection Code
prefix = "linux"
if platform == "linux" or platform == "linux32":
	prefix = "linux"
elif platform == "darwin":
	prefix = "osx"
elif platform == "win32":
	prefix = "win64"

dllpath = path.join(mz.__path__[0], prefix + '_' + 'libmozart.so')
lib = CDLL(dllpath)

def getMatrix1D(degree):
	"""
	Get FEM matrices on the reference domain

	Paramters
		- ``degree`` : degree of polynomial

	Returns
		- ``M_R`` : 
		- ``S_R`` : 
		- ``D_R`` :  

	"""
	M_R = None
	S_R = None
	D_R = None
	if degree is 1:
		M_R = np.array([[ 2,  1], [ 1, 2]]) / 3.
		S_R = np.array([[ 1, -1], [-1, 1]]) / 2.
		D_R = np.array([[-1,  1], [-1, 1]]) / 2.
	if degree is 2:
		M_R = np.array([[ 4,  2, -1], [ 2, 16,  2], [-1,  2, 4]]) / 15.
		S_R = np.array([[ 7, -8,  1], [-8, 16, -8], [ 1, -8, 7]]) / 6.
		D_R = np.array([[-3,  4, -1], [-1,  0,  1], [ 1, -4, 3]]) / 2.

	return (M_R, S_R, D_R)

def nJacobiP(x=np.array([0]), alpha=0, beta=0, degree=0):
	"""
	Evaluate normalized Jacobi polynomial of type alpha, beta > -1 at point x for order n
	(the special case of alpha = beta = 0, knows as the normalized Legendre polynomial)

	Paramters
		- ``x`` (``float64 array``) : variable x
		- ``alpha`` (``int32``) : superscript alpha of normalized Jacobi polynomial
		- ``beta`` (``int32``) : superscript beta of normalized Jacobi polynomial
		- ``degree`` (``int32``) : Polynomial degree

	Returns
		- ``P`` (``float64 array``) : the value of degree-th order normalized Jacobi polynomial at x
	
	Example
		>>> N = 2
		>>> x = np.array([-1, 0, 1])
		>>> from mozart.poisson.solve import nJacobiP
		>>> P = nJacobiP(x,0,0,N)
		>>> print(P)
		array([ 1.58113883, -0.79056942,  1.58113883])	
	"""
	Pn = np.zeros((degree+1,x.size),float)
	Pn[0,:] = np.sqrt(2**(-alpha-beta-1) * np.math.gamma(alpha+beta+2) / ((np.math.gamma(alpha + 1) * np.math.gamma(beta + 1))))

	if degree == 0:
		P = Pn
	else:
		Pn[1,:] = np.multiply(Pn[0,:]*np.sqrt((alpha+beta+3.0)/((alpha+1)*(beta+1))),((alpha+beta+2)*x+(alpha-beta)))/2
		a_n = 2.0/(2+alpha+beta)*np.sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
		for n in range(2,degree+1):
			anew=2.0/(2*n+alpha+beta)*np.sqrt(n*(n+alpha+beta)*(n+alpha)*(n+beta)/((2*n+alpha+beta-1)*(2*n+alpha+beta+1)))
			b_n=-(alpha**2-beta**2)/((2*(n-1)+alpha+beta)*(2*(n-1)+alpha+beta+2.0))
			Pn[n,:]=(np.multiply((x-b_n),Pn[n-1,:])-a_n*Pn[n-2,:])/anew;
			a_n=anew;

	P = Pn[degree,:]
	return P

def DnJacobiP(x=np.array([0]), alpha=0, beta=0, degree=0):
	"""
	Evaluate the derivative of the normalized Jacobi polynomial of type alpha, beta > -1 at point x for order n

	Paramters
		- ``x`` (``float64 array``) : variable x
		- ``alpha`` (``int32``) : superscript alpha of normalized Jacobi polynomial
		- ``beta`` (``int32``) : superscript beta of normalized Jacobi polynomial
		- ``degree`` (``int32``) : Polynomial degree

	Returns
		- ``dP`` (``float64 array``) : the value of the derivative of the normalized Jacobi polynomial at x
									   according to alpha, beta, degree
	
	Example
		>>> N = 2
		>>> x = np.array([-1, 0, 1])
		>>> from mozart.poisson.solve import DnJacobiP
		>>> dP = DnJacobiP(x,0,0,N)
		>>> print(dP)
		array([-4.74341649,  0.        ,  4.74341649])	
	"""
	length = x.size
	dP = np.zeros(length,float)
	if degree == 0:
		dP[:] = 0;
	else:
		dP[:] = np.sqrt(degree*(degree+alpha+beta+1))*nJacobiP(x,alpha+1,beta+1,degree-1);
	return dP

def one_dim(c4n, n4e, n4Db, f, u_D, degree = 1):
	"""
	Computes the coordinates of nodes and elements.
	
	Parameters
		- ``c4n`` (``float64 array``) : coordinates
		- ``n4e`` (``int32 array``) : nodes for elements
{'fraction_non_alphanumeric': 0.21848225214198286, 'fraction_numerical': 0.05538555691554468, 'mean_word_length': 3.6100141043723553, 'pattern_counts': {'":': 6, '<': 3, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3064912', 'n_tokens_mistral': 3070, 'n_tokens_neox': 2768, 'n_words': 688}