Skip to content
Snippets Groups Projects
Commit 4aa0c7f4 authored by Anthony Royer's avatar Anthony Royer
Browse files

Remove demos and missing contrib/bessel

parent b875f17b
No related branches found
No related tags found
No related merge requests found
Pipeline #8132 passed
Stage: build
Stage: test
Show whitespace changes
Inline Side-by-side
Showing
with 8198 additions and 1071 deletions
contrib/bessel/Bessel.cpp 0 → 100644
+ 272
− 0
View file @ 4aa0c7f4
// GetDP - Copyright (C) 1997-2020 P. Dular and C. Geuzaine, University of Liege
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include "Bessel.h"
#include <iostream>
#if defined(HAVE_NO_FORTRAN)
static void zbesj_(double*, double*, double*, int*, int*, double*,
double*, int*, int*)
{
std::cerr << "Bessel functions require Fortran compiler";
}
static void zbesk_(double*, double*, double*, int*, int*, double*,
double*, int*, int*)
{
std::cerr << "Bessel functions require Fortran compiler";
}
static void zbesy_(double*, double*, double*, int*, int*, double*,
double*, int*, double*, double*, int*)
{
std::cerr << "Bessel functions require Fortran compiler";
}
static void zbesh_(double*, double*, double*, int*, int*, int*,
double*, double*, int*, int*)
{
std::cerr << "Bessel functions require Fortran compiler";
}
static void zairy_(double*, double*, int*, int*, double*,
double*, int*, int*)
{
std::cerr << "Bessel functions require Fortran compiler";
}
#else
#if defined(HAVE_UNDERSCORE)
#define zbesj_ zbesj
#define zbesk_ zbesk
#define zbesy_ zbesy
#define zbesh_ zbesh
#define zairy_ zairy
#endif
extern "C" {
void zbesj_(double*, double*, double*, int*, int*, double*,
double*, int*, int*);
void zbesk_(double*, double*, double*, int*, int*, double*,
double*, int*, int*);
void zbesy_(double*, double*, double*, int*, int*, double*,
double*, int*, double*,
double*, int*);
void zbesh_(double*, double*, double*, int*, int*, int*, double*,
double*, int*, int*);
void zairy_(double*, double*, int*, int*, double*,
double*, int*, int*);
}
#endif
static int BesselError(int ierr, const char *str)
{
static int warn=0;
switch(ierr){
case 0 :
return 0;
case 1 :
std::cerr << "Input error in " << str;
return BESSEL_ERROR_INPUT;
case 2 :
return BESSEL_OVERFLOW;
case 3 :
if(!warn){
std::cout << "Half machine accuracy lost in " << str << " (large argument or order)";
warn = 1;
}
return BESSEL_HALF_ACCURACY;
case 4 :
std::cerr << "Complete loss of significance in " << str << " (argument or order too large)";
return BESSEL_NO_ACCURACY;
case 5 :
std::cerr << "Failed to converge in " << str;
return BESSEL_NO_CONVERGENCE;
default:
std::cout << "Unknown Bessel status in " << str << " (" << ierr << ")";
return ierr;
}
}
// First kind Bessel functions
int BesselJn(double n, int num, double x, double *val)
{
int nz = 0, ierr = 0, kode = 1;
double xi = 0.0;
double* ji = new double[num];
zbesj_(&x, &xi, &n, &kode, &num, val, ji, &nz, &ierr) ;
delete[] ji;
return BesselError(ierr, "BesselJn");
}
int BesselJnComplex(double n, int num, double xr, double xi, double *valr, double *vali)
{
int nz = 0, ierr = 0, kode = 1;
zbesj_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr) ;
return BesselError(ierr, "BesselJnComplex");
}
int BesselKnComplex(double n, int num, double xr, double xi, double *valr, double *vali)
{
int nz = 0, ierr = 0, kode = 1;
zbesk_(&xr, &xi, &n, &kode, &num, valr, vali, &nz, &ierr) ;
return BesselError(ierr, "BesselKnComplex");
}
int BesselSphericalJn(double n, int num, double x, double *val)
{
int ierr = BesselJn(n+0.5, num, x, val);
double coef = sqrt(0.5*M_PI/x);
for(int i = 0; i < num; i++){
val[i] *= coef;
}
return BesselError(ierr, "BesselSphericalJn");
}
int BesselAltSphericalJn(double n, int num, double x, double *val)
{
int ierr = BesselJn(n+0.5, num, x, val);
double coef = sqrt(0.5*M_PI*x);
for(int i = 0; i < num; i++){
val[i] *= coef;
}
return BesselError(ierr, "BesselAltSphericalJn");
}
// Second kind Bessel functions
int BesselYn(double n, int num, double x, double *val)
{
int nz = 0, ierr = 0, kode = 1;
double xi = 0.0;
double* yi = new double[num];
double* auxyr = new double[num];
double* auxyi = new double[num];
zbesy_(&x, &xi, &n, &kode, &num, val, yi, &nz, auxyr, auxyi, &ierr);
delete[] yi;
delete[] auxyr;
delete[] auxyi;
return BesselError(ierr, "BesselYn");
}
int BesselSphericalYn(double n, int num, double x, double *val)
{
int ierr = BesselYn(n+0.5, num, x, val);
double coef = sqrt(0.5*M_PI/x);
for(int i = 0; i < num; i++){
val[i] *= coef;
}
return BesselError(ierr, "BesselSphericalYn");
}
int BesselAltSphericalYn(double n, int num, double x, double *val)
{
int ierr = BesselYn(n+0.5, num, x, val);
double coef = sqrt(0.5*M_PI*x);
for(int i = 0; i < num; i++){
val[i] *= coef;
}
return BesselError(ierr, "BesselAltSphericalYn");
}
// Hankel functions (type = 1 or 2)
int BesselHn(int type, double n, int num, double x, std::complex<double> *val)
{
int nz = 0, ierr = 0, kode = 1;
double* hr = new double[num];
double* hi = new double[num];
double xi = 0.0;
zbesh_(&x, &xi, &n, &kode, &type, &num, hr, hi, &nz, &ierr);
for(int i=0; i < num; i++){
val[i] = std::complex<double>(hr[i], hi[i]);
}
delete[] hr;
delete[] hi;
return BesselError(ierr, "BesselHn");
}
int BesselSphericalHn(int type, double n, int num, double x, std::complex<double> *val)
{
int ierr = BesselHn(type, n+0.5, num, x, val);
double coef = sqrt(0.5*M_PI/x);
for(int i = 0; i < num; i++){
val[i] *= coef;
}
return BesselError(ierr, "BesselSphericalHn");
}
int BesselAltSphericalHn(int type, double n, int num, double x, std::complex<double> *val)
{
int ierr = BesselHn(type, n+0.5, num, x, val);
double coef = sqrt(0.5*M_PI*x);
for(int i = 0; i < num; i++){
val[i] *= coef;
}
return BesselError(ierr, "BesselAltSphericalHn");
}
// Airy Function and derivative (id=1)
int AiryComplex(double xr, double xi, int id, double *valr, double *vali)
{
int nz = 0, ierr = 0, kode = 1;
zairy_(&xr, &xi, &id, &kode, valr, vali, &nz, &ierr);
//val = std::complex<double>(valr, vali);
return BesselError(ierr, "AiryComplex");
}
// Utilities for backward compatibility
double Spherical_j_n(int n, double x)
{
double res;
BesselSphericalJn(n, 1, x, &res);
return res;
}
double AltSpherical_j_n(int n, double x)
{
double res;
BesselAltSphericalJn(n, 1, x, &res);
return res;
}
double Spherical_y_n(int n, double x)
{
double res;
BesselSphericalYn(n, 1, x, &res);
return res;
}
double AltSpherical_y_n(int n, double x)
{
double res;
BesselAltSphericalYn(n, 1, x, &res);
return res;
}
contrib/bessel/Bessel.h 0 → 100644
+ 45
− 0
View file @ 4aa0c7f4
// GetDP - Copyright (C) 1997-2020 P. Dular and C. Geuzaine, University of Liege
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues.
#ifndef BESSEL_H
#define BESSEL_H
#include <cmath>
#include <complex>
#define BESSEL_ERROR_INPUT 1
#define BESSEL_OVERFLOW 2
#define BESSEL_HALF_ACCURACY 3
#define BESSEL_NO_ACCURACY 4
#define BESSEL_NO_CONVERGENCE 5
// These routines provide a C++ interface to the Fortran Bessel
// functions from Donald E. Amos (Sandia National Laboratories)
int BesselJn(double n, int num, double x, double *val);
int BesselSphericalJn(double n, int num, double x, double *val);
int BesselAltSphericalJn(double n, int num, double x, double *val);
int BesselJnComplex(double n, int num, double xr, double xi, double *valr, double *vali);
int BesselKnComplex(double n, int num, double xr, double xi, double *valr, double *vali);
int BesselYn(double n, int num, double x, double *val);
int BesselSphericalYn(double n, int num, double x, double *val);
int BesselAltSphericalYn(double n, int num, double x, double *val);
int BesselHn(int type, double n, int num, double x, std::complex<double> *val);
int BesselSphericalHn(int type, double n, int num, double x, std::complex<double> *val);
int BesselAltSphericalHn(int type, double n, int num, double x, std::complex<double> *val);
int AiryComplex(double xr, double xi, int id, double *valr, double *vali);
// Utilities for backward compatibility
double Spherical_j_n(int n, double x);
double AltSpherical_j_n(int n, double x);
double Spherical_y_n(int n, double x);
double AltSpherical_y_n(int n, double x);
#endif
contrib/bessel/BesselLib.f 0 → 100644
+ 7866
− 0
View file @ 4aa0c7f4
Source diff could not be displayed: it is too large. Options to address this: view the blob.
contrib/bessel/CMakeLists.txt 0 → 100644
+ 15
− 0
View file @ 4aa0c7f4
# GetDP - Copyright (C) 1997-2020 P. Dular and C. Geuzaine, University of Liege
#
# See the LICENSE.txt file for license information. Please report all
# issues on https://gitlab.onelab.info/getdp/getdp/issues.
set(SRC
Bessel.cpp
)
if(ENABLE_BESSEL)
list(APPEND SRC BesselLib.f)
endif(ENABLE_BESSEL)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
append_gmshfem_contrib(Numeric "${SRC};${HDR}")
examples/helmholtzFlow/Duct_ABCs/CMakeLists.txt deleted 100644 → 0
+ 0
− 8
View file @ b875f17b
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
project(demo CXX)
include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake)
add_executable(demo main.cpp ${EXTRA_INCS})
target_link_libraries(demo ${EXTRA_LIBS})
examples/helmholtzFlow/Duct_ABCs/main.cpp deleted 100644 → 0
+ 0
− 267
View file @ b875f17b
#include "GmshFem.h"
#include "Formulation.h"
#include "gmsh.h"
#include "AnalyticalFunction.h"
using namespace gmshfem;
using namespace gmshfem::common;
using namespace gmshfem::problem;
using namespace gmshfem::domain;
using namespace gmshfem::field;
using namespace gmshfem::function;
using namespace gmshfem::analytics;
using namespace gmshfem::post;
using namespace gmshfem::equation;
void meshDuct(const double L, const double H, const double h, const int MeshOrder)
{
msg::info << " - Duct of length L = " << L << msg::endl;
msg::info << " - Duct of heigth H = " << H << msg::endl;
gmsh::model::add("geometry");
// duct definition
gmsh::model::geo::addPoint(L, H, 0., h, 6);
gmsh::model::geo::addPoint(0, H, 0., h, 7);
gmsh::model::geo::addPoint(0, 0, 0., h, 8);
gmsh::model::geo::addPoint(L, 0, 0., h, 9);
gmsh::model::geo::addLine(6, 7, 5);
gmsh::model::geo::addLine(7, 8, 6);
gmsh::model::geo::addLine(8, 9, 7);
gmsh::model::geo::addLine(9, 6, 8);
gmsh::model::geo::addCurveLoop({5,6,7,8}, 1);
gmsh::model::geo::addPlaneSurface({1}, 1);
gmsh::model::geo::mesh::setTransfiniteSurface(1);
gmsh::model::geo::mesh::setRecombine(1,1);
// physical groups
gmsh::model::addPhysicalGroup(2, {1}, 1);
gmsh::model::setPhysicalName(2, 3, "omega");
gmsh::model::addPhysicalGroup(1, {5}, 2);
gmsh::model::setPhysicalName(1, 2, "gammaTop");
gmsh::model::addPhysicalGroup(1, {6}, 3);
gmsh::model::setPhysicalName(1, 3, "gammaLeft");
gmsh::model::addPhysicalGroup(1, {7}, 4);
gmsh::model::setPhysicalName(1, 4, "gammaBottom");
gmsh::model::addPhysicalGroup(1, {8}, 5);
gmsh::model::setPhysicalName(1, 5, "gammaRight");
gmsh::model::addPhysicalGroup(0, {6}, 6);
gmsh::model::setPhysicalName(0, 6, "gammaTopCorner");
gmsh::model::addPhysicalGroup(0, {9}, 7);
gmsh::model::setPhysicalName(0, 7, "gammaBottomCorner");
// generate mesh
gmsh::model::geo::synchronize();
gmsh::model::mesh::generate(2);
gmsh::model::mesh::setOrder(MeshOrder);
gmsh::option::setNumber("Mesh.ElementOrder",MeshOrder);
gmsh::model::mesh::setOrder(MeshOrder);
}
int main(int argc, char **argv)
{
//*****
// Problem declaration
//*****
GmshFem gmshFem(argc, argv);
const double pi = 3.14159265358979323846264338327950288;
const std::complex< double > im(0., 1.);
// numerical parameters
int FEMorder = 8;
gmshFem.userDefinedParameter(FEMorder, "FEMorder");
std::string gauss = "Gauss"+std::to_string(2*FEMorder+1);
gmshFem.userDefinedParameter(gauss, "gauss");
// meshing
double L = 0.5; // duct length
gmshFem.userDefinedParameter(L, "L");
double H = 0.25; // duct height
gmshFem.userDefinedParameter(H, "H");
double h = 1./40; // meshsize
gmshFem.userDefinedParameter(h, "h");
meshDuct(L, H, h, 1);
// duct problem specifications
std::string problem = "hard";
gmshFem.userDefinedParameter(problem, "problem");
double w = 30; // frequency
gmshFem.userDefinedParameter(w, "omega");
std::vector< unsigned int > n = {3}; // single mode
// std::vector< unsigned int > n = {0,2,4,6,8}; // multiple modes modes
std::vector< double > A(n.size(),1.); // select amplitudes
int N_modes = n.size();
msg::info << " There are " << N_modes << " input modes " << msg::endl;
// Choose ABC type
std::string abcName = "Pade";
gmshFem.userDefinedParameter(abcName, "abcName");
// mean flow
const double M = 0.8;
double beta = sqrt(1-M*M);
std::vector< std::complex< double > > KX(N_modes);
std::vector< double > KY(N_modes), SIGMA(N_modes);
for(unsigned int j = 0; j < n.size(); ++j) {
KY[j] = n[j]*pi/H;
SIGMA[j] = w*w-beta*beta*KY[j]*KY[j];
if ( SIGMA[j] >=0 ) {
KX[j] = (1/(beta*beta))*(-M*w + sqrt(SIGMA[j]));
msg::info << " Mode " << n[j] << " is propagative " << msg::endl;
if ( w>(beta*KY[j]) && w<KY[j] ) {
msg::info << " Mode " << n[j] << " is inverse upstream ! " << msg::endl;
}
}
else {
KX[j] = (1/(beta*beta))*(-M*w - im* sqrt(abs(SIGMA[j])));
msg::info << " Mode " << n[j] << " is evanescent " << msg::endl;
}
msg::info << " - propagating wavenumber, kx = " << KX[j] << " of mode " << n[j] << msg::endl;
}
float PtsByWl = 2*pi*FEMorder*(1+M) / w / h;
msg::info << "********************************" << msg::endl;
msg::info << " - convected Helmholtz duct problem " << msg::endl;
msg::info << " - wave number = " << w << msg::endl;
msg::info << " - FEM basis order = " << FEMorder << "" << msg::endl;
msg::info << " - Approximate dofs by wavelength = " << PtsByWl << "" << msg::endl;
msg::info << "********************************" << msg::endl;
// declare formulation and FEM domains
std::vector< FieldInterface< std::complex< double > > * > fieldBucket;
Formulation< std::complex< double > > formulation("helmholtzflow");
Domain omega(2, 1);
Domain gammaTop(1, 2);
Domain gammaLeft(1, 3);
Domain gammaBottom(1, 4);
Domain gammaRight(1, 5);
Domain gammaTopCorner(0, 6);
Domain gammaBottomCorner(0, 7);
Field< std::complex< double >, Form::Form0 > v("v", omega | gammaLeft | gammaRight | gammaTop | gammaBottom, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
// define analytical solution
Function< std::complex< double >, Degree::Degree0 > *solution = nullptr;
if(problem == "soft") {
v.addConstraint(gammaTop, 0.);
v.addConstraint(gammaBottom, 0.);
solution = new AnalyticalFunction< helmholtz2D::DuctModeSolutionMulti< std::complex< double > > >(w, M, H, A, n, 0., 0., 1);
}
else {
solution = new AnalyticalFunction< helmholtz2D::DuctModeSolutionMulti< std::complex< double > > >(w, M, H, A, n, 0., 0., 0);
}
solution->activateMemorization();
formulation.galerkin(vector< std::complex< double > >(beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss);
formulation.galerkin(vector< std::complex< double > >(0., 1., 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omega, gauss);
formulation.galerkin(vector< std::complex< double > >(im*w*M, 0., 0.) * grad(dof(v)), tf(v), omega, gauss);
formulation.galerkin(-im*w*M*dof(v), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss);
formulation.galerkin(- w * w * dof(v), tf(v), omega, gauss);
// input BC
ScalarFunction< std::complex< double > > Input_modal_sum=0.;
for(unsigned int j = 0; j < n.size(); ++j) {
if(problem == "soft"){
Input_modal_sum = Input_modal_sum + A[j] * KX[j] * sin( KY[j] * y< std::complex< double > >() );
}
else if(problem == "hard"){
Input_modal_sum = Input_modal_sum + A[j] * KX[j] * cos( KY[j] * y< std::complex< double > >() );
}
else {
msg::error << " this problem is not defined " << msg::endl;
}
}
formulation.galerkin(-beta*beta * im * Input_modal_sum , tf(v), gammaLeft, gauss);
formulation.galerkin(-im * w * M * dof(v), tf(v), gammaLeft, gauss);
// output BC
if(abcName == "ABC-0") {
msg::info << "Use Sommerfeld ABC." << msg::endl;
const double angle = 0.;
std::complex< double > wM = w*(M-exp(im*angle/2.)) / (beta*beta); //angle=0 -> -k/(1+M);
formulation.galerkin(-beta*beta * im * wM * dof(v), tf(v), gammaRight, gauss);
formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss);
}
else if(abcName == "ABC-2") {
msg::info << "Use second order ABC." << msg::endl;
const double angle = 0.;
// contribution to the mass and rigidity matrices
std::complex< double > wM = w*(M-cos(angle/2.)) / (beta*beta);
std::complex< double > wK = -exp(-im*angle/2.)/(2*w);
formulation.galerkin(-beta*beta * im * wM * dof(v), tf(v), gammaRight, gauss);
formulation.galerkin( beta*beta * im * wK * grad(dof(v)), grad(tf(v)), gammaRight, gauss);
formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss);
/* use DtN eigenvalues instead of Laplace Beltrami operator
std::complex< double > kT = -k + KY[0]*KY[0]/(2*k);
formulation.galerkin(-beta*beta * im * kT * dof(v), tf(v), gammaRight, gauss);
*/
}
else if(abcName == "Pade") {
int padeOrder = 4;
gmshFem.userDefinedParameter(padeOrder, "padeOrder");
double angle = -pi/4.;
gmshFem.userDefinedParameter(angle, "angle");
msg::info << "Use Pade ABC of order " << padeOrder << " with angle " << angle << " rad"<< msg::endl;
const double Np = 2. * padeOrder + 1.;
const std::complex< double > exp1 = std::complex<double>(std::cos(angle),std::sin(angle));
const std::complex< double > exp2 = std::complex<double>(std::cos(angle/2.),std::sin(angle/2.));
const std::complex< double > coef = 2./Np;
std::vector< std::complex< double > > c(padeOrder, 0.);
for(int i = 0; i < padeOrder; ++i) {
c[i] = std::tan((i + 1) * pi / Np);
c[i] *= c[i];
}
// define the auxiliary fields
std::vector< Field< std::complex< double >, Form::Form0 >* > phi;
Field< std::complex< double >, Form::Form0 >* phis;
for(int i = 0; i < padeOrder; ++i) {
phis = new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i), gammaRight | gammaTopCorner | gammaBottomCorner, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
if(problem == "soft") { // impose Dirichlet constraints on the auxiliary fields
phis->addConstraint(gammaTopCorner,0.);
phis->addConstraint(gammaBottomCorner,0.);
}
phi.push_back(phis);
fieldBucket.push_back(phi.back());
}
// write the augmented weak form
formulation.galerkin(im * w * exp2 * dof(v), tf(v), gammaRight, gauss);
for(int i = 0; i < padeOrder; ++i) {
// boundary integral terms relating the auxiliary fields
formulation.galerkin(im * w * exp2 * coef * c[i] * dof(*phi[i]), tf(v), gammaRight, gauss);
formulation.galerkin(im * w * exp2 * coef * c[i] * dof(v), tf(v), gammaRight, gauss);
// coupling of the auxiliary equations
formulation.galerkin( grad(dof(*phi[i])), grad(tf(*phi[i])), gammaRight, gauss);
formulation.galerkin(- (w*w)/(beta*beta) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), gammaRight, gauss);
formulation.galerkin(- (w*w)/(beta*beta) * exp1 * (c[i] + 1.) * dof(v), tf(*phi[i]), gammaRight, gauss);
}
}
else if(abcName == "DtN") {
msg::info << "Use analytical DtN." << msg::endl;
formulation.galerkin(beta*beta * im * KX[0] * dof(v), tf(v), gammaRight, gauss);
formulation.galerkin(im * w * M * dof(v), tf(v), gammaRight, gauss);
}
formulation.pre();
formulation.assemble();
formulation.solve();
save(+v, omega, "v");
save(*solution, omega, "v_exact");
save(*solution - v, omega, "error");
std::complex< double > num = integrate(pow(abs(*solution - v), 2), omega, gauss);
std::complex< double > den = integrate(pow(abs(*solution), 2), omega, gauss);
msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl;
for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
delete fieldBucket[i];
}
return 0;
}
examples/helmholtzFlow/Duct_PMLs/CMakeLists.txt deleted 100644 → 0
+ 0
− 8
View file @ b875f17b
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
project(demo CXX)
include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake)
add_executable(demo main.cpp ${EXTRA_INCS})
target_link_libraries(demo ${EXTRA_LIBS})
examples/helmholtzFlow/Duct_PMLs/main.cpp deleted 100644 → 0
+ 0
− 216
View file @ b875f17b
#include "GmshFem.h"
#include "Formulation.h"
#include "gmsh.h"
#include "AnalyticalFunction.h"
using namespace gmshfem;
using namespace gmshfem::common;
using namespace gmshfem::problem;
using namespace gmshfem::domain;
using namespace gmshfem::field;
using namespace gmshfem::function;
using namespace gmshfem::analytics;
using namespace gmshfem::post;
using namespace gmshfem::equation;
int main(int argc, char **argv)
{
GmshFem gmshFem(argc, argv);
double L = 0.5; // duct length
gmshFem.userDefinedParameter(L, "L");
double H = 0.25; // duct height
gmshFem.userDefinedParameter(H, "H");
int Npml = 4; // number of PML layers
gmshFem.userDefinedParameter(Npml, "Npml");
double h = 2*sqrt(L*H)/50; // meshsize
gmshFem.userDefinedParameter(h, "h");
int MeshOrder = 1;
gmshFem.userDefinedParameter(MeshOrder, "MeshOrder");
gmsh::model::add("helmholtzflow");
// physical duct
gmsh::model::geo::addPoint(L, H, 0., h, 6);
gmsh::model::geo::addPoint(0, H, 0., h, 7);
gmsh::model::geo::addPoint(0, 0, 0., h, 8);
gmsh::model::geo::addPoint(L, 0, 0., h, 9);
gmsh::model::geo::addLine(6, 7, 5);
gmsh::model::geo::addLine(7, 8, 6);
gmsh::model::geo::addLine(8, 9, 7);
gmsh::model::geo::addLine(9, 6, 8);
gmsh::model::geo::addCurveLoop({5,6,7,8}, 1);
gmsh::model::geo::addPlaneSurface({1}, 1);
gmsh::model::geo::mesh::setTransfiniteSurface(1);
gmsh::model::geo::mesh::setRecombine(1,1);
// pml extension
gmsh::model::geo::addPoint(L+Npml*h, H, 0., h, 10);
gmsh::model::geo::addPoint(L+Npml*h, 0, 0., h, 11);
gmsh::model::geo::addLine(6, 10, 12);
gmsh::model::geo::addLine(10, 11, 13);
gmsh::model::geo::addLine(11, 9, 14);
gmsh::model::geo::addCurveLoop({12,13,14,8}, 2);
gmsh::model::geo::addPlaneSurface({2}, 2);
gmsh::model::geo::mesh::setTransfiniteSurface(2);
gmsh::model::geo::mesh::setRecombine(2,2);
// physical groups
gmsh::model::addPhysicalGroup(2, {1}, 1);
gmsh::model::setPhysicalName(2, 3, "omega");
gmsh::model::addPhysicalGroup(1, {5}, 2);
gmsh::model::setPhysicalName(1, 2, "gammaTop");
gmsh::model::addPhysicalGroup(1, {6}, 3);
gmsh::model::setPhysicalName(1, 3, "gammaLeft");
gmsh::model::addPhysicalGroup(1, {7}, 4);
gmsh::model::setPhysicalName(1, 4, "gammaBottom");
gmsh::model::addPhysicalGroup(1, {8}, 5);
gmsh::model::setPhysicalName(1, 5, "gammaRight");
// pml groups
gmsh::model::addPhysicalGroup(2, {2}, 6);
gmsh::model::setPhysicalName(2, 6, "omegaPml");
gmsh::model::addPhysicalGroup(1, {12}, 7);
gmsh::model::setPhysicalName(1, 7, "gammaTopPml");
gmsh::model::addPhysicalGroup(1, {13}, 8);
gmsh::model::setPhysicalName(1, 8, "gammaRightPml");
gmsh::model::addPhysicalGroup(1, {14}, 9);
gmsh::model::setPhysicalName(1, 9, "gammaBottomPml");
// generate mesh
gmsh::model::geo::synchronize();
gmsh::model::mesh::generate(2);
//gmsh::write("m.msh");
gmsh::model::mesh::setOrder(MeshOrder);
//*****
// Problem declaration
//*****
int FEMorder = 4;
gmshFem.userDefinedParameter(FEMorder, "FEMorder");
std::string gauss = "Gauss10";
gmshFem.userDefinedParameter(gauss, "gauss");
std::string problem = "hard";
gmshFem.userDefinedParameter(problem, "problem");
std::string PmlType = "Alternative"; // Lorentz, Classic, Alternative
gmshFem.userDefinedParameter(PmlType, "PmlType");
double pi = 3.14159265359;
const double k = 70;
const double M = 0.8;
const int n = 3;
const std::complex< double > im(0., 1.);
double ky = n*pi/H;
double beta = sqrt(1-M*M);
double sigma = k*k-beta*beta*ky*ky;
std::complex< double > kx;
if (sigma >= 0) {
kx = (1/(beta*beta))*(-M*k + sqrt(sigma));
msg::info << " propagative mode " << msg::endl;
}
else {
kx = (1/(beta*beta))*(-M*k - im* sqrt(abs(sigma)));
msg::info << " evanescent mode " << msg::endl;
}
msg::info << "********************************" << msg::endl;
msg::info << " - convected Helmholtz duct problem " << msg::endl;
msg::info << " - wave number = " << k << msg::endl;
msg::info << " - propagating wavenumber = " << kx << msg::endl;
if ( k>(beta*ky) && k<ky )
{
msg::info << " - inverse upstream mode ! " << msg::endl;
}
unsigned int pointsByWl = (2*pi*FEMorder/k) / h;
msg::info << " - FEM basis order = " << FEMorder << "" << msg::endl;
msg::info << " - Approximate dofs by wavelength = " << pointsByWl << "" << msg::endl;
msg::info << "********************************" << msg::endl;
Formulation< std::complex< double > > formulation("helmholtzflow");
Domain omega(2, 1);
Domain omegaPml(2, 6);
Domain gammaTop(1, 2);
Domain gammaLeft(1, 3);
Domain gammaBottom(1, 4);
Domain gammaRight(1, 5);
Domain gammaTopPml(1, 7);
Domain gammaBottomPml(1, 9);
Field< std::complex< double >, Form::Form0 > v("v", omega | omegaPml | gammaLeft | gammaRight | gammaTop | gammaBottom | gammaTopPml | gammaBottomPml, FunctionSpaceTypeForm0::HierarchicalH1, FEMorder);
// define analytical solution
Function< std::complex< double >, Degree::Degree0 > *solution = nullptr;
if(problem == "soft"){
v.addConstraint(gammaTop, 0.);
v.addConstraint(gammaBottom, 0.);
v.addConstraint(gammaBottomPml, 0.);
v.addConstraint(gammaTopPml, 0.);
solution = new AnalyticalFunction< helmholtz2D::DuctModeSolution< std::complex< double > > >(k, M, H, 0., 0., n, 1);
}
else {
solution = new AnalyticalFunction< helmholtz2D::DuctModeSolution< std::complex< double > > >(k, M, H, 0., 0., n, 0);
}
solution->activateMemorization();
formulation.galerkin(vector< std::complex< double > >(beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss);
formulation.galerkin(vector< std::complex< double > >(0., 1., 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omega, gauss);
formulation.galerkin(vector< std::complex< double > >(im*k*M, 0., 0.) * grad(dof(v)), tf(v), omega, gauss);
formulation.galerkin(-im*k*M*dof(v), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omega, gauss);
formulation.galerkin(- k * k * dof(v), tf(v), omega, gauss);
// input BC
if(problem == "soft"){
formulation.galerkin(-beta*beta * im * kx * (sin(ky * y< std::complex< double > >())) , tf(v), gammaLeft, gauss);
}
else{
formulation.galerkin(-beta*beta * im * kx * (cos(ky * y< std::complex< double > >())) , tf(v), gammaLeft, gauss);
}
formulation.galerkin(-im * k * M * dof(v), tf(v), gammaLeft, gauss);
msg::info << "Use a PML" << msg::endl;
msg::info << "PML formulation: " << PmlType << msg::endl;
const double Lpml = L+Npml*h;
double Sigma0 = 4*beta*beta;
// Bermudez function
ScalarFunction< std::complex< double > > SigmaX = Sigma0/(Lpml - x< std::complex< double > >() );
ScalarFunction< std::complex< double > > gammaX = 1-(im/k)*SigmaX;
ScalarFunction< std::complex< double > > gammaXinv = 1. / gammaX;
// Different types of PML formulations
if (PmlType == "Classic"){
formulation.galerkin(vector< std::complex< double > >(gammaXinv*beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(0., gammaX, 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(im*k*M, 0., 0.) * grad(dof(v)), tf(v), omegaPml, gauss);
formulation.galerkin(-im*k*M*dof(v), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(- gammaX*k*k * dof(v), tf(v), omegaPml, gauss);
}
else if(PmlType == "Lorentz"){
formulation.galerkin(vector< std::complex< double > >(gammaXinv*beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(0., gammaX, 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(gammaXinv*im*k*M, 0., 0.) * grad(dof(v)), tf(v), omegaPml, gauss);
formulation.galerkin(-gammaXinv*im*k*M*dof(v), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(- k*k*(gammaX-gammaXinv*M*M)/(beta*beta) * dof(v), tf(v), omegaPml, gauss);
}
else if(PmlType == "Alternative"){
formulation.galerkin(vector< std::complex< double > >(gammaXinv*beta*beta, 0., 0.) * grad(dof(v)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(0., gammaX, 0.) * grad(dof(v)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(v)), omegaPml, gauss);
formulation.galerkin(- gammaX*k*k/(beta*beta) * dof(v), tf(v), omegaPml, gauss);
}
formulation.pre();
formulation.assemble();
formulation.solve();
save(+v, omega, "v");
save(+v, omegaPml, "vpml");
save(*solution, omega, "v_exact");
save(*solution - v, omega, "error");
std::complex< double > num = integrate(pow(abs(*solution - v), 2), omega, gauss);
std::complex< double > den = integrate(pow(abs(*solution), 2), omega, gauss);
msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl;
return 0;
}
examples/helmholtzFlow/FreeField/CMakeLists.txt deleted 100644 → 0
+ 0
− 8
View file @ b875f17b
cmake_minimum_required(VERSION 3.0 FATAL_ERROR)
project(demo CXX)
include(${CMAKE_CURRENT_SOURCE_DIR}/../../../demos.cmake)
add_executable(demo main.cpp ${EXTRA_INCS})
target_link_libraries(demo ${EXTRA_LIBS})
examples/helmholtzFlow/FreeField/main.cpp deleted 100644 → 0
+ 0
− 564
View file @ b875f17b
#include "GmshFem.h"
#include "Formulation.h"
#include "gmsh.h"
#include "AnalyticalFunction.h"
using namespace gmshfem;
using namespace gmshfem::common;
using namespace gmshfem::problem;
using namespace gmshfem::domain;
using namespace gmshfem::field;
using namespace gmshfem::function;
using namespace gmshfem::analytics;
using namespace gmshfem::post;
using namespace gmshfem::equation;
struct Edge
{
Domain gamma;
Domain corner[2];
};
struct HelmholtzDomain
{
Domain omega;
Domain gammaPhy;
Domain gammaPml;
Domain omegaPml;
Domain source;
Domain corners;
std::vector< Domain > corner;
std::vector< Edge > edge;
};
//***************************************
// GEOMETRY
//***************************************
void meshCircle(GmshFem &gmshFem, HelmholtzDomain &domains, const double R, const double Rpml, const double xs, const double ys, const double lc, bool withPml)
{
gmsh::model::add("geometry");
gmsh::model::geo::addPoint(0., 0., 0., lc, 1); // Center
gmsh::model::geo::addPoint(xs, ys, 0., lc, 100); // Source point
// physical circle
gmsh::model::geo::addPoint(R, 0., 0., lc, 2);
gmsh::model::geo::addPoint(0, R, 0., lc, 3);
gmsh::model::geo::addPoint(-R, 0., 0., lc, 4);
gmsh::model::geo::addPoint(0, -R, 0., lc, 5);
gmsh::model::geo::addCircleArc(2, 1, 3, 1);
gmsh::model::geo::addCircleArc(3, 1, 4, 2);
gmsh::model::geo::addCircleArc(4, 1, 5, 3);
gmsh::model::geo::addCircleArc(5, 1, 2, 4);
gmsh::model::geo::addCurveLoop({1,2,3,4}, 1);
// surface
gmsh::model::geo::addPlaneSurface({1}, 1);
if(withPml) {
gmsh::model::geo::addPoint(Rpml, 0., 0., lc, 6);
gmsh::model::geo::addPoint(0, Rpml, 0., lc, 7);
gmsh::model::geo::addPoint(-Rpml, 0., 0., lc, 8);
gmsh::model::geo::addPoint(0, -Rpml, 0., lc, 9);
gmsh::model::geo::addCircleArc(6, 1, 7, 5);
gmsh::model::geo::addCircleArc(7, 1, 8, 6);
gmsh::model::geo::addCircleArc(8, 1, 9, 7);
gmsh::model::geo::addCircleArc(9, 1, 6, 8);
gmsh::model::geo::addCurveLoop({5,6,7,8}, 2);
gmsh::model::geo::addPlaneSurface({1,2}, 2);
}
gmsh::model::geo::synchronize();
gmsh::model::mesh::embed(0,{100},2,1);
// physicals
gmsh::model::addPhysicalGroup(2, {1}, 1);
gmsh::model::setPhysicalName(2, 1, "omega");
gmsh::model::addPhysicalGroup(1, {1,2,3,4}, 3);
gmsh::model::setPhysicalName(1, 3, "gammaPhy");
gmsh::model::addPhysicalGroup(0, {100}, 100);
gmsh::model::setPhysicalName(0, 100, "source");
if(withPml) {
gmsh::model::addPhysicalGroup(2, {2}, 2);
gmsh::model::setPhysicalName(2, 2, "omegaPml");
gmsh::model::addPhysicalGroup(1, {5,6,7,8}, 4);
gmsh::model::setPhysicalName(1, 4, "gammaPml");
}
// generate mesh
gmsh::model::geo::synchronize();
gmsh::model::mesh::generate(2);
// gmsh::model::mesh::setOrder(2); // T6 elem
// gmsh::write("m.msh");
// Domains
domains.omega = Domain(2, 1);
domains.gammaPhy = Domain(1, 3);
domains.source = Domain(0, 100);
if(withPml) {
domains.omegaPml = Domain(2, 2);
domains.gammaPml = Domain(1, 4);
}
}
void meshSquare(GmshFem &gmshFem, HelmholtzDomain &domains, const double L, const double Lpml, const double xs, const double ys, const double lc, bool withPml)
{
gmsh::model::add("geometry");
gmsh::model::geo::addPoint(0., 0., 0., lc, 1); // Center
gmsh::model::geo::addPoint(xs, ys, 0., lc, 100); // Source point
// polygone
gmsh::model::geo::addPoint(L, L, 0., lc, 2);
gmsh::model::geo::addPoint(-L, L, 0., lc, 3);
gmsh::model::geo::addPoint(-L, -L, 0., lc, 4);
gmsh::model::geo::addPoint(L, -L, 0., lc, 5);
gmsh::model::geo::addLine(2, 3, 1);
gmsh::model::geo::addLine(3, 4, 2);
gmsh::model::geo::addLine(4, 5, 3);
gmsh::model::geo::addLine(5, 2, 4);
gmsh::model::geo::addCurveLoop({1,2,3,4}, 1);
gmsh::model::geo::addPlaneSurface({1}, 1);
if(withPml) {
gmsh::model::geo::addPoint(Lpml, Lpml, 0., lc, 6);
gmsh::model::geo::addPoint(-Lpml, Lpml, 0., lc, 7);
gmsh::model::geo::addPoint(-Lpml, -Lpml, 0., lc, 8);
gmsh::model::geo::addPoint(Lpml, -Lpml, 0., lc, 9);
gmsh::model::geo::addLine(6, 7, 5);
gmsh::model::geo::addLine(7, 8, 6);
gmsh::model::geo::addLine(8, 9, 7);
gmsh::model::geo::addLine(9, 6, 8);
gmsh::model::geo::addCurveLoop({5,6,7,8}, 2);
gmsh::model::geo::addPlaneSurface({1,2}, 2);
}
gmsh::model::geo::synchronize();
gmsh::model::mesh::embed(0,{100},2,1);
// physicals
gmsh::model::addPhysicalGroup(2, {1}, 1);
gmsh::model::setPhysicalName(2, 1, "omega");
gmsh::model::addPhysicalGroup(1, {1,2,3,4}, 3);
gmsh::model::setPhysicalName(1, 3, "gammaPhy");
gmsh::model::addPhysicalGroup(0, {100}, 100);
gmsh::model::setPhysicalName(0, 100, "source");
// Physical edges
gmsh::model::addPhysicalGroup(1, {1}, 20);
gmsh::model::setPhysicalName(1, 20, "gammaPhy_top");
gmsh::model::addPhysicalGroup(1, {2}, 30);
gmsh::model::setPhysicalName(1, 30, "gammaPhy_left");
gmsh::model::addPhysicalGroup(1, {3}, 40);
gmsh::model::setPhysicalName(1, 40, "gammaPhy_bottom");
gmsh::model::addPhysicalGroup(1, {4}, 50);
gmsh::model::setPhysicalName(1, 50, "gammaPhy_right");
// Physical domain corners
gmsh::model::addPhysicalGroup(0, {2}, 2);
gmsh::model::setPhysicalName(0, 2, "UpRightCorner");
gmsh::model::addPhysicalGroup(0, {3}, 3);
gmsh::model::setPhysicalName(0, 3, "UpLeftCorner");
gmsh::model::addPhysicalGroup(0, {4}, 4);
gmsh::model::setPhysicalName(0, 4, "BottomLeftCorner");
gmsh::model::addPhysicalGroup(0, {5}, 5);
gmsh::model::setPhysicalName(0, 5, "BottomRightCorner");
if(withPml) {
gmsh::model::addPhysicalGroup(2, {2}, 2);
gmsh::model::setPhysicalName(2, 2, "omegaPml");
gmsh::model::addPhysicalGroup(1, {5,6,7,8}, 4);
gmsh::model::setPhysicalName(1, 4, "gammaPml");
}
// generate mesh
gmsh::model::geo::synchronize();
gmsh::model::mesh::generate(2);
//gmsh::model::mesh::setOrder(2); // T6 elem
// gmsh::write("m.msh");
// Domains
domains.omega = Domain(2, 1);
domains.gammaPhy = Domain(1, 3);
domains.source = Domain(0, 100);
if(withPml) {
domains.omegaPml = Domain(2, 2);
domains.gammaPml = Domain(1, 4);
}
domains.edge.resize(4);
domains.corner.resize(4);
for(int i = 0; i < 4; ++i) {
domains.edge[i].gamma = Domain(1, 10*(i + 2) );
domains.edge[i].corner[0] = Domain(0, i + 2 );
domains.edge[i].corner[1] = Domain(0, (i + 1) % 4 + 2 );
domains.corner[i] = Domain(0, i + 2);
domains.corners |= domains.edge[i].corner[0];
domains.gammaPhy |= domains.edge[i].gamma;
}
}
//***************************************
// FORMULATION
//***************************************
int main(int argc, char **argv)
{
GmshFem gmshFem(argc, argv);
// geometry parameters
int geometry = 0; // 0-circle, 1-square
gmshFem.userDefinedParameter(geometry, "geometry");
double R = 1.0; // Circle radius or square of size [-R,R] x [-R,R]
gmshFem.userDefinedParameter(R, "R");
// Point source location
double xs = 0.3;
double ys = 0.2;
gmshFem.userDefinedParameter(xs, "xs");
gmshFem.userDefinedParameter(ys, "ys");
msg::info << "Point source located at xs=" << xs << ", ys=" << ys << msg::endl;
// physical parameters
const double pi = 3.14159265358979323846264338327950288;
const double k = 3 * pi; // free field wavenumber
const std::complex< double > im(0., 1.);
double M = 0.; // Mach number
gmshFem.userDefinedParameter(M, "M");
double theta = pi/4; // mean flow orientation
gmshFem.userDefinedParameter(theta, "theta");
msg::info << "Mach number " << M << " oriented at theta=" << theta << " rad" << msg::endl;
// numerical parameters
int order = 4; // FEM shape function order
gmshFem.userDefinedParameter(order, "order");
std::string gauss = "Gauss10"; // increase Nr of Gauss points when order > 6
double lc = 0.03; // mesh size
gmshFem.userDefinedParameter(lc, "lc");
// Choose ABC - exterior boundary condition
std::string abcName = "Pade"; // available choices: ABC-0, ABC-2, Pade, pml
gmshFem.userDefinedParameter(abcName, "abcName");
int Npml;
double Rpml;
if (abcName == "pml") { // PML parameters
Npml = 4; // number of layers
gmshFem.userDefinedParameter(Npml, "Npml");
Rpml = R+Npml*lc; // extended domain
}
HelmholtzDomain domains;
switch (geometry) {
case 0:
meshCircle(gmshFem, domains, R, Rpml, xs, ys, lc, abcName == "pml");
break;
case 1:
meshSquare(gmshFem, domains, R, Rpml, xs, ys, lc, abcName == "pml");
break;
default:
msg::error << "geometry not available ! " << msg::endl;
exit(0);
break;
}
float pointsByWl = (2*pi*order*(1-abs(M)))/(lc*k);
msg::info << " - FEM basis order = " << order << "" << msg::endl;
msg::info << " - Smallest number of dofs by wavelength = " << pointsByWl << "" << msg::endl;
if (pointsByWl < 6) {
msg::warning << " - Less than 6 points per wavelength ! " << msg::endl;
}
// Mach number projections
double Mx = M*std::cos(theta);
double My = M*std::sin(theta);
// Mach-velocity vector
VectorFunction< std::complex< double > > MM = vector< std::complex< double > > (Mx,My,0.);
// Normal and tangential projections
ScalarFunction< std::complex< double > > Mn = vector< std::complex< double > > (Mx,My,0.)*normal< std::complex< double > >();
ScalarFunction< std::complex< double > > Mt = vector< std::complex< double > > (Mx,My,0.)*tangent< std::complex< double > >();
// Parameters for Inverse Lorentz transformation
double beta = sqrt(1-M*M); // Jacobian
double Alphax = 1 + Mx*Mx/(beta*(1+beta));
double Alphay = 1 + My*My/(beta*(1+beta));
double K = Mx*My/(beta*(1+beta));
// Tensor of the Inverse Lorentz transformation
TensorFunction< std::complex< double > > Linv = tensor< std::complex <double > > (beta*Alphay, -beta*K, 0., -beta*K, beta*Alphax,0.,0.,0.,0.);
std::vector< FieldInterface< std::complex< double > > * > fieldBucket;
Formulation< std::complex< double > > formulation("helmholtzflow");
Field< std::complex< double >, Form::Form0 > u("u", domains.omega | domains.omegaPml | domains.gammaPhy | domains.gammaPml | domains.corners | domains.source, FunctionSpaceTypeForm0::HierarchicalH1, order);
// Define analytic solution
Function< std::complex< double >, Degree::Degree0 > *solution = nullptr;
solution = new AnalyticalFunction< helmholtz2D::MonopoleFreeField< std::complex< double > > >(k, M, theta, 1., xs, ys, 0., 0.);
solution->activateMemorization();
// convected Helmholz weak form
formulation.galerkin(vector< std::complex< double > >(1-Mx*Mx, -Mx*My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), domains.omega, gauss);
formulation.galerkin(vector< std::complex< double > >(-Mx*My, 1-My*My, 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omega, gauss);
formulation.galerkin(- k * k * dof(u), tf(u), domains.omega, gauss);
formulation.galerkin(dof(u), vector< std::complex< double > >(-im*k*Mx,-im*k*My, 0.) * grad(tf(u)), domains.omega, gauss);
formulation.galerkin(vector< std::complex< double > >(im*k*Mx, im*k*My, 0.) * grad(dof(u)), tf(u), domains.omega, gauss);
// Source point
formulation.galerkin(-1.,tf(u),domains.source,gauss);
// Exact DtN: beta_n**2 * kx = - Mn * (k0 - Mt*ky) + sqrt( k0**2 -2*k0*Mt*ky - (1-abs(M))^2 ky**2 )
if(abcName == "ABC-0") { // zeroth order Taylor approximation
// boundary contributions
formulation.galerkin(im*k*Mn * dof(u), tf(u), domains.gammaPhy, gauss);
formulation.galerkin(Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
// ABC
msg::info << "Use zeroth order ABC" << msg::endl;
formulation.galerkin(im*k*(1-Mn) * dof(u), tf(u), domains.gammaPhy, gauss);
formulation.galerkin(-Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
}
else if (abcName == "ABC-2") {
// boundary contributions
formulation.galerkin(im*k*Mn * dof(u), tf(u), domains.gammaPhy, gauss);
formulation.galerkin(Mn*Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
// ABC
msg::info << "Use second order ABC" << msg::endl;
// zeroth order contribution
formulation.galerkin(im*k*(1-Mn) * dof(u), tf(u), domains.gammaPhy, gauss);
// first order contribution
formulation.galerkin( (1-Mn) * Mt * tangent< std::complex< double > >() * grad(dof(u)), dof(u), domains.gammaPhy, gauss);
// second order contribution
formulation.galerkin(- im * (beta*beta) /(2*k) * grad(dof(u)), grad(tf(u)), domains.gammaPhy, gauss);
}
else if (abcName == "Pade") {
int padeOrder = 2;
gmshFem.userDefinedParameter(padeOrder, "padeOrder");
double angle = -pi/4.;
gmshFem.userDefinedParameter(angle, "angle");
msg::info << "Use Pade ABC of order " << padeOrder << " with angle " << angle << " rad"<< msg::endl;
const double Np = 2. * padeOrder + 1.;
const std::complex< double > exp1 = std::complex<double>(std::cos(angle),std::sin(angle));
const std::complex< double > exp2 = std::complex<double>(std::cos(angle/2.),std::sin(angle/2.));
const std::complex< double > coef = 2./Np;
std::vector< std::complex< double > > c(padeOrder, 0.);
for(int i = 0; i < padeOrder; ++i) {
c[i] = std::tan((i + 1) * pi / Np);
c[i] *= c[i];
}
if(geometry == 0) {
// define the auxiliary fields
std::vector< Field< std::complex< double >, Form::Form0 >* > phi;
for(int i = 0; i < padeOrder; ++i) {
phi.push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i), domains.gammaPhy, FunctionSpaceTypeForm0::HierarchicalH1, order));
fieldBucket.push_back(phi.back());
}
// write the augmented weak form - approximation of the square-root
formulation.galerkin( im * k * exp2 * dof(u), tf(u), domains.gammaPhy, gauss);
for(int i = 0; i < padeOrder; ++i) {
// boundary integral terms relating the auxiliary fields
formulation.galerkin( im * k * exp2 * coef * c[i] * dof(*phi[i]), tf(u), domains.gammaPhy, gauss);
formulation.galerkin( im * k * exp2 * coef * c[i] * dof(u), tf(u), domains.gammaPhy, gauss);
// coupling of the auxiliary equations
formulation.galerkin(- (beta*beta) * grad(dof(*phi[i])), grad(tf(*phi[i])), domains.gammaPhy, gauss);
formulation.galerkin(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[i])), tf(*phi[i]), domains.gammaPhy, gauss);
formulation.galerkin( (k*k) * (exp1 * c[i] + 1.) * dof(*phi[i]), tf(*phi[i]), domains.gammaPhy, gauss);
formulation.galerkin( (k*k) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[i]), domains.gammaPhy, gauss);
}
}
else if(geometry == 1) {
// define the auxilary fields
std::vector< Field< std::complex< double >, Form::Form0 > * > phi[4];
for(int x = 0; x < 4; ++x) {
for(int i = 0; i < padeOrder; ++i) {
phi[x].push_back(new Field< std::complex< double >, Form::Form0 >("phi_" + std::to_string(i) + "^" + std::to_string(x), domains.edge[x].gamma | domains.edge[x].corner[0] | domains.edge[x].corner[1], FunctionSpaceTypeForm0::HierarchicalH1,order));
fieldBucket.push_back(phi[x].back());
}
}
formulation.galerkin( im * k * exp2 * dof(u), tf(u), domains.gammaPhy, gauss);
for(int x = 0; x < 4; ++x) {
for(int i = 0; i < padeOrder; ++i) {
// Boundary conditions
formulation.galerkin( im * k * exp2 * coef * c[i] * dof(*phi[x][i]), tf(u), domains.edge[x].gamma, gauss);
formulation.galerkin( im * k * exp2 * coef * c[i] * dof(u), tf(u), domains.edge[x].gamma, gauss);
// Auxiliary equations
formulation.galerkin(- (beta*beta) * grad(dof(*phi[x][i])), grad(tf(*phi[x][i])), domains.edge[x].gamma, gauss);
formulation.galerkin(- 2. * im * k * Mt * tangent< std::complex< double > >() * grad(dof(*phi[x][i])), tf(*phi[x][i]), domains.edge[x].gamma, gauss);
formulation.galerkin( (k*k) * (exp1 * c[i] + 1.) * dof(*phi[x][i]), tf(*phi[x][i]), domains.edge[x].gamma, gauss);
formulation.galerkin( (k*k) * exp1 * (c[i] + 1.) * dof(u), tf(*phi[x][i]), domains.edge[x].gamma, gauss);
}
}
bool withCorner = true;
gmshFem.userDefinedParameter(withCorner, "withCorner");
std::string CornerStrategy = "Sommerfeld"; // Sommerfeld or HABC
gmshFem.userDefinedParameter(CornerStrategy, "CornerStrategy");
if(withCorner) {
msg::info << "with corner treatment" << msg::endl;
msg::info << "use " << CornerStrategy << " condition at corners" << msg::endl;
// define the corner fields
std::vector< std::vector< std::vector< Field< std::complex< double >, Form::Form0 > * > > > corner_phi(4);
for(int k = 0; k < 4; ++k) {
corner_phi[k].resize(padeOrder);
for(int i = 0; i < padeOrder; ++i) {
for(int j = 0; j < padeOrder; ++j) {
corner_phi[k][i].push_back(new Field< std::complex< double >, Form::Form0 >("corner_phi_(" + std::to_string(i) + ", " + std::to_string(j) + ")^(" + std::to_string(k) + ")", domains.corner[k], FunctionSpaceTypeForm0::HierarchicalH1,order));
fieldBucket.push_back(corner_phi[k][i].back());
}
}
}
for(int ci = 0; ci < 4; ++ci) {
int x = ci;
int y = (ci - 1 < 0 ? 4 - 1 : ci - 1);
for(int i = 0; i < padeOrder; ++i) {
if (CornerStrategy == "Sommerfeld"){
// strategy 1 - Sommerfeld condition at corners
formulation.galerkin(- im * k * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss);
formulation.galerkin(- im * k * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
}
else if (CornerStrategy == "HABC") {
// strategy 2 - HABC at corners
formulation.galerkin(- im * k * exp2 * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss);
formulation.galerkin(- im * k * exp2 * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
for(int j = 0; j < padeOrder; ++j) {
// Auxiliary equations for corner
formulation.galerkin( (exp1*c[i] + exp1*c[j] + 1.) * dof(*corner_phi[ci][i][j]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss);
formulation.galerkin( exp1 * (c[j] + 1.) * dof(*phi[x][i]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss);
formulation.galerkin( exp1 * (c[i] + 1.) * dof(*phi[y][j]), tf(*corner_phi[ci][i][j]), domains.corner[ci], gauss);
// Corner conditions
formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*phi[x][i]), tf(*phi[x][i]), domains.corner[ci], gauss);
formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*phi[y][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*corner_phi[ci][i][j]), tf(*phi[x][i]), domains.corner[ci], gauss);
formulation.galerkin(- im * k * exp2 * coef * c[j] * dof(*corner_phi[ci][j][i]), tf(*phi[y][i]), domains.corner[ci], gauss);
}
}
else {
msg::error << "Corner strategy not available !" << msg::endl;
exit(0);
}
}
}
}
}
}
else if (abcName == "pml") {
msg::info << "Use a PML with " << Npml << " layers" << msg::endl;
double Sigma0 = beta;
ScalarFunction< std::complex< double > > det_J;
TensorFunction< std::complex< double > > J_PML_inv_T;
const double Wpml = Rpml - R;
if (geometry == 0) {
msg::info << "stabilized PML - circular domain" << msg::endl;
msg::info << "use unbounded absorbing function with sigma_0 = " << Sigma0 << msg::endl;
ScalarFunction< std::complex< double > > cosT = x< std::complex< double > >() / r2d< std::complex< double > >();
ScalarFunction< std::complex< double > > sinT = y< std::complex< double > >() / r2d< std::complex< double > >();
ScalarFunction< std::complex< double > > dampingProfileR = Sigma0 / (Wpml - (r2d< std::complex< double > >() - R));
ScalarFunction< std::complex< double > > dampingProfileInt = -Sigma0*ln((Wpml - (r2d< std::complex< double > >() - R)) / Wpml);
ScalarFunction< std::complex< double > > gamma = 1. - im * dampingProfileR / k;
ScalarFunction< std::complex< double > > gamma_hat = 1. - im * (1. / r2d< std::complex< double > >()) * dampingProfileInt / k;
det_J = gamma * gamma_hat;
J_PML_inv_T = tensor< std::complex <double > >(cosT/gamma,sinT/gamma,0.,-sinT/gamma_hat,cosT/gamma_hat,0.,0.,0.,0.);
}
else if (geometry == 1) {
msg::info << "stabilized PML - square domain" << msg::endl;
msg::info << "use unbounded absorbing function with sigma_0 = " << Sigma0 << msg::endl;
ScalarFunction< std::complex< double > > SigmaX = heaviside(abs(x< std::complex< double > >()) - R ) * Sigma0/(Rpml - abs(x< std::complex< double > >()) );
ScalarFunction< std::complex< double > > gammaX = 1-(im/k)*SigmaX;
ScalarFunction< std::complex< double > > SigmaY = heaviside(abs(y< std::complex< double > >()) - R ) * Sigma0/(Rpml - abs(y< std::complex< double > >()) );
ScalarFunction< std::complex< double > > gammaY = 1-(im/k)*SigmaY;
det_J = gammaX * gammaY;
J_PML_inv_T = tensor< std::complex <double > >(1./gammaX,0.,0.,0.,1./gammaY,0.,0.,0.,0.);
/*
msg::info << "classical PML - square domain - unstable ! " << msg::endl;
formulation.galerkin(vector< std::complex< double > >((gammaY/gammaX)*(1-Mx*Mx), -Mx*My, 0.) * grad(dof(u)), vector< std::complex< double > >(1., 0., 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(-Mx*My, (gammaX/gammaY)*(1-My*My), 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin(dof(u), vector< std::complex< double > >(-gammaY*im*k*Mx,-gammaX*im*k*My, 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin(vector< std::complex< double > >(gammaY*im*k*Mx, gammaX*im*k*My, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss);
formulation.galerkin(- k*k* gammaX*gammaY * dof(u), tf(u), domains.omegaPml, gauss);
*/
/*
msg::info << "Alternative stable PML - square domain only " << msg::endl;
ScalarFunction< std::complex< double > > Ax = (gammaY/gammaX)*(1-Mx*Mx);
std::complex< double > Cx =-im*k*Mx*( 1 - (My*My)/(1-Mx*Mx) ) / (beta*beta);
double Ka = -Mx*My/(1-Mx*Mx);
ScalarFunction< std::complex< double > > Ay = (gammaX/gammaY)*( 1 - (My*My)/(1-Mx*Mx) );
std::complex< double > Cy = -im*k*My/(beta*beta);
// modified bilinear term in dx
formulation.galerkin( Ax * vector< std::complex< double > >(1. , Ka, 0.) * grad(dof(u)), vector< std::complex< double > >(1., Ka, 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin( Ax * Cx * dof(u), vector< std::complex< double > >(1., Ka, 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin( Ax * (-Cx) * vector< std::complex< double > >(1. , Ka, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss);
formulation.galerkin( Ax * Cx * (-Cx) *dof(u), tf(u), domains.omegaPml, gauss);
// modified bilinear term in dy
formulation.galerkin( Ay * vector< std::complex< double > >(0. , 1., 0.) * grad(dof(u)), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin( Ay * Cy * dof(u), vector< std::complex< double > >(0., 1., 0.) * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin( Ay * vector< std::complex< double > >(0. , -Cy, 0.) * grad(dof(u)), tf(u), domains.omegaPml, gauss);
formulation.galerkin( Ay * Cy * (-Cy) *dof(u), tf(u), domains.omegaPml, gauss);
//
formulation.galerkin(- k*k/(beta*beta) * gammaX*gammaY * dof(u), tf(u), domains.omegaPml, gauss);
*/
}
// build vector and matrix for the weak form
VectorFunction< std::complex< double > > J_PML_inv_T_M = J_PML_inv_T * MM;
TensorFunction< std::complex< double > > J_PML_Linv = J_PML_inv_T * Linv;
// stabilized PML weak form - valid for a general domain
formulation.galerkin( det_J * J_PML_Linv*grad(dof(u)) , J_PML_Linv*grad(tf(u)) , domains.omegaPml, gauss);
formulation.galerkin(+ k*k/(beta*beta) * det_J * J_PML_inv_T_M * dof(u) , J_PML_inv_T_M * tf(u), domains.omegaPml, gauss);
formulation.galerkin(+ im*k/beta* det_J * J_PML_Linv * grad(dof(u)), J_PML_inv_T_M * tf(u), domains.omegaPml, gauss);
formulation.galerkin(- im*k/beta* det_J * J_PML_inv_T_M * dof(u), J_PML_Linv * grad(tf(u)), domains.omegaPml, gauss);
formulation.galerkin(- k*k/(beta*beta) * det_J * dof(u), tf(u), domains.omegaPml, gauss);
}
else {
msg::error << "ABC Type not available !" << msg::endl;
exit(0);
}
// solve
formulation.pre();
formulation.assemble();
formulation.solve();
// compute the projection of the analytic solution on the FEM basis - best FEM approximation
Formulation< std::complex< double > > projection("helmholtzflow");
Field< std::complex< double >, Form::Form0 > uP("uP", domains.omega, FunctionSpaceTypeForm0::HierarchicalH1, order);
projection.galerkin( -dof(uP), tf(uP), domains.omega, gauss);
projection.galerkin( *solution, tf(uP), domains.omega, gauss);
msg::info << "Compute best approximation... " << msg::endl;
projection.pre();
projection.assemble();
projection.solve();
// save results and compute errors
save(+u, domains.omega, "u");
save(+uP, domains.omega, "uP");
if (abcName == "pml"){
save(+u, domains.omegaPml, "u_pml");
}
save(*solution, domains.omega, "u_exact");
save(*solution - u, domains.omega, "error");
std::complex< double > num = integrate(pow(abs(*solution - u), 2), domains.omega, gauss);
std::complex< double > numP = integrate(pow(abs(*solution - uP), 2), domains.omega, gauss);
std::complex< double > den = integrate(pow(abs(*solution), 2), domains.omega, gauss);
msg::info << "L_2 error = " << 100.*sqrt(num / den) << " %" << msg::endl;
msg::info << "Best L_2 error = " << 100.*sqrt(numP / den) << " %" << msg::endl;
msg::warning << "L_2 error should exclude neighbourhood elements to the point source" << msg::endl;
for(unsigned int i = 0; i < fieldBucket.size(); ++i) {
delete fieldBucket[i];
}
return 0;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment