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
[email protected]
https://mail.gna.org/listinfo/getfem-users
