-
Notifications
You must be signed in to change notification settings - Fork 47
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Draft: add basis transformation procedure #984
base: main
Are you sure you want to change the base?
Changes from 1 commit
70aa74e
c7b1451
74b8495
1eebdd9
6ade65d
4c3e40c
233a8d6
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -1219,6 +1219,58 @@ int CeedBasisDestroy(CeedBasis *basis) { | |
return CEED_ERROR_SUCCESS; | ||
} | ||
|
||
/** | ||
@brief Hale and Trefethen strip transformation | ||
|
||
@param[in] Q Number of quadrature points | ||
@param[in] rho sum of semiminor and major axis | ||
@param[out] g conformal map | ||
@param[out] g_prime derivative of conformal map g | ||
|
||
@return An error code: 0 - success, otherwise - failure | ||
|
||
@ref Utility | ||
**/ | ||
int CeedHaleTrefethenStripMap(CeedInt Q, CeedInt rho, CeedScalar *g, CeedScalar *g_prime) { | ||
CeedScalar tau, d, C, PI2, PI = 4.0*atan(1.0); | ||
tau = PI / log(rho); | ||
d = .5 + 1. / (exp(tau * PI) + 1.0) | ||
PI2 = PI / 2.0; | ||
for (int i = 0; i < Q; i++) { | ||
u = (PI * i) / Q; | ||
// Unscaled function of u | ||
g_u = log(1.0 + exp(-tau * (PI2 + u))) - log(1.0 + exp(-tau * (PI2 - u))) + d * tau * u; | ||
g_prime_u = 1.0 / (exp(tau * (PI2 + u)) + 1) + 1.0 / (exp(tau * (PI2 - u)) + 1) - d; | ||
// Normalizing factor and scaled function | ||
C = 1.0 / g_u; | ||
g = C * g_u; | ||
g_prime = -tau * C ; | ||
} | ||
return CEED_ERROR_SUCCESS; | ||
} | ||
|
||
/** | ||
@brief Transform quadrature by applying a smooth mapping = (g, g_prime) | ||
|
||
@param Q Number of quadarture points | ||
|
||
@return An error code: 0 - success, otherwise - failure | ||
|
||
@ref Utility | ||
**/ | ||
int CeedTransformQuadrature(CeeddInt Q, CeedScalar *q_weight_1d) { | ||
CeedScalar points, m_points, m_weights; | ||
if (!mapping) { | ||
if (!q_weight_1d) | ||
points = Q; | ||
else | ||
points = Q; | ||
weights = q_weight_1d; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Note that the indentation here does not show the scope. Your editor should be able to re-indent to make the actual scope of the |
||
} | ||
|
||
return CEED_ERROR_SUCCESS; | ||
} | ||
|
||
/** | ||
@brief Construct a Gauss-Legendre quadrature | ||
|
||
|
@@ -1270,6 +1322,8 @@ int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, | |
q_ref_1d[i] = -xi; | ||
q_ref_1d[Q-1-i]= xi; | ||
} | ||
// Call transformed Gauss-Legendre quadrature | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. You probably can't just call it here because you can't change the behavior of the existing interfaces, but you'll add a new interface that creates a Gauss quadrature and then transforms it. |
||
|
||
return CEED_ERROR_SUCCESS; | ||
} | ||
|
||
|
@@ -1342,6 +1396,8 @@ int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, | |
q_ref_1d[i] = -xi; | ||
q_ref_1d[Q-1-i]= xi; | ||
} | ||
// Call transformed Gauss-Legendre-Lobatto quadrature | ||
|
||
return CEED_ERROR_SUCCESS; | ||
} | ||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If this is a public function, it'll need to be declared in a public header.