Re: [Getfem-users] Add a new FEM with polynomial and global functions

Dear Zaim Yassine

It would be useful to know which version of GetFEM++ you are currently using. If you are using version 5.1 or later I can help you at least regarding the use of global basis functions. The standard way I am using them is by deriving a class from getfem::global_function_simple, lets call it global_function_custom, where I redefine the following methods

virtual scalar_type val(const base_node &pt) const
virtual void grad(const base_node &pt, base_small_vector &g) const
virtual void hess(const base_node &pt, base_matrix &h) const
virtual bool is_in_support(const base_node &pt) const
virtual void bounding_box(base_node &bmin, base_node &bmax) const

where the coordinates of the point pt are global coordinates.

Then I save all instances of my derived class, in a vector of shared pointers

std::vector<getfem::pglobal_function> basis_funcs;
for (...)
basis_funcs.push_back(std::make_shared<global_function_custom>(...));

Then I define my mesh_fem as

getfem::mesh_fem_global_function mf(m,1);
mf.set_functions(basis_funcs, mim);

and I am ready to use mf.

The difference between global fem and standard fem is that the basis functions in the latter case are defined in terms of the local/element coordinates and there is an additional geometric transformation between global and local coordinates.

I hope this helps you a bit.

Best regards

Kostas

On Sun, Nov 20, 2016 at 6:38 PM, Yassine ZAIM <address@hidden> wrote:

Hello dear GetFEM users,
I am trying to add a new FEM with some global functions in addition to some polynomial functions. I found this discussion about how I can do it:
https://mail.gna.org/public/getfem-users/2014-04/msg00034.html#content
However when I look in the file "getfem_mesh_fem_global_function.cc/h" to inspire. I don't found where the basis functions are defined? as the functions:
std::stringstream s
for (int i = 0; i<size_basis; ++i) p->base()[i] =bgeot::read_base_poly(dim,s);
for the polynomial case; And also where the name of FEM is defined? as the function: add_suffix("Name", fem_element); in the polynomial case.

I tried to add my element by programming a similar class of
template <class FUNC> class fem : public virtual_fem {...};
in which I have determined explicitly the value of base_value, grad_base_value and hess_base_value. I inherited from this class to define my element in the file getfem_fem.cc. I define the basis and DOF by the functions:
//############# code to add the basis and DOF ###############//
base_[i] = polynomial or global basis function;//like base_value in annex
add_node(DOF, Point);// corresponding to each basis function

I know that for this method I could program a class with the functions eval() and derivative(). But in my case I defined the functions base_value, grad_base_value and hess_base_value explicitly without need of these methods (I think). You can see my class in the annex. By this way I get a bad result.

I hope that I was clear, and I will be thankful for your help of how I can add my element correctly.

///########## Annex ##############///
class MyFUNC : public virtual_fem {
protected :
std::vector<opt_long_scalar_type> base_;

public :
/// Gives the array of basic functions (components).
const std::vector<opt_long_scalar_type> &base(void) const { return base_; }
std::vector<opt_long_scalar_type> &base(void) { return base_; }

/** Evaluates at point x, all base functions and returns the result in
t(nb_base,target_dim) */
void base_value(const base_node &z, base_tensor &t) const {
//scalar_type res = 0;
bgeot::multi_index mi(2);
mi[1] = target_dim(); mi[0] = short_type(nb_base(0));
t.adjust_sizes(mi);
base_tensor::iterator it = t.begin();
scalar_type x = *z.begin();//z[0];
scalar_type y = *z.end();
*it = bgeot::to_scalar(x*y); ++it;
*it = bgeot::to_scalar((1-x)*y); ++it;
*it = bgeot::to_scalar((1-x)*(1-y)); ++it;
*it = bgeot::to_scalar(x*(1-y)); ++it;
*it = bgeot::to_scalar((32/1281)*sqrt(2)*(1-x)*sqrt(1-x)*(12*x-66*x*x+(143/2)*x*x*x)); ++it;
*it = bgeot::to_scalar((32/1281)*sqrt(2)*(1-y)*sqrt(1-y)*(12*y+66*y*y+(143/2)*y*y*y));
}
/** Evaluates at point x, the gradient of all base functions w.r.t. the
reference element directions 0,..,dim-1 and returns the result in
t(nb_base,target_dim,dim) */
void grad_base_value(const base_node &z, base_tensor &t) const {
bgeot::multi_index mi(3);
dim_type n = dim();
mi[2] = n; mi[1] = target_dim(); mi[0] = short_type(nb_base(0));
t.adjust_sizes(mi);
base_tensor::iterator it = t.begin();
scalar_type x = *z.begin();
scalar_type y = *z.end();
*it = bgeot::to_scalar(y); ++it;
*it = bgeot::to_scalar(-y); ++it;
*it = bgeot::to_scalar(-(1-y) ); ++it;
*it = bgeot::to_scalar((1-y)); ++it;
*it = bgeot::to_scalar((32/1281)*sqrt(2)*(-3/2*sqrt(1-x)*(12*x-66*x*x+(143/2)*x*x*x)+(1-x)*sqrt(1-x)*(12-132*x+(429/2)*x*x))); ++it;
*it = bgeot::to_scalar(0); ++it;

*it = bgeot::to_scalar(x); ++it;
*it = bgeot::to_scalar((1-x)); ++it;
*it = bgeot::to_scalar(-(1-x)); ++it;
*it = bgeot::to_scalar(-x); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar((32/1281)*sqrt(2)*(-3/2*sqrt(1-y)*(12*y-66*y*y+(143/2)*y*y*y)+(1-y)*sqrt(1-y)*(12-132*y+(429/2)*y*y)));
}
/** Evaluates at point x, the hessian of all base functions w.r.t. the
reference element directions 0,..,dim-1 and returns the result in
t(nb_base,target_dim,dim,dim) */
void hess_base_value(const base_node &z, base_tensor &t) const {
bgeot::multi_index mi(4);
dim_type n = dim();
mi[3] = n; mi[2] = n; mi[1] = target_dim();
mi[0] = short_type(nb_base(0));
t.adjust_sizes(mi);
base_tensor::iterator it = t.begin();
scalar_type x = *z.begin();
scalar_type y = *z.end();

*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar((32/1281)*sqrt(2)*(3/4*(1/sqrt(1-x))*(12*x-66*x*x+(143/2)*x*x*x)-3*sqrt(1-x)*(12-132*x+(429/2)*x*x)+(1-x)*sqrt(1-x)*(-132+429*x))); ++it;
*it = bgeot::to_scalar(0); ++it;

*it = bgeot::to_scalar(1); ++it;
*it = bgeot::to_scalar(-1); ++it;
*it = bgeot::to_scalar(1); ++it;
*it = bgeot::to_scalar(-1); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;

*it = bgeot::to_scalar(1); ++it;
*it = bgeot::to_scalar(-1); ++it;
*it = bgeot::to_scalar(1); ++it;
*it = bgeot::to_scalar(-1); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;

*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar(0); ++it;
*it = bgeot::to_scalar((32/1281)*sqrt(2)*(3/4*(1/sqrt(1-y))*(12*y-66*y*y+(143/2)*y*y*y)-3*sqrt(1-y)*(12-132*y+(429/2)*y*y)+(1-y)*sqrt(1-y)*(-132+429*y)));
}
};

--
ZAIM Yassine
PhD Student in Applied Mathematics

_______________________________________________
Getfem-users mailing list
address@hidden
https://mail.gna.org/listinfo/getfem-users

From:	Konstantinos Poulios
Subject:	Re: [Getfem-users] Add a new FEM with polynomial and global functions
Date:	Sun, 20 Nov 2016 22:51:35 +0100