25 template <
typename Derived>
26 void secDerToHessian(
const gsEigen::DenseBase<Derived> & secDers,
28 gsMatrix<typename Derived::Scalar> & hessian)
30 const index_t sz = dim*(dim+1)/2;
31 const gsAsConstMatrix<typename Derived::Scalar>
32 ders(secDers.derived().data(), sz, secDers.size() / sz );
33 hessian.resize(dim*dim, ders.cols() );
38 hessian = secDers.transpose();
41 hessian.row(0)=ders.row(0);
43 hessian.row(2)=ders.row(2);
44 hessian.row(3)=ders.row(1);
47 hessian.row(0)=ders.row(0);
49 hessian.row(1)=ders.row(3);
51 hessian.row(2)=ders.row(4);
52 hessian.row(4)=ders.row(1);
54 hessian.row(5)=ders.row(5);
55 hessian.row(8)=ders.row(2);
66 fs(&_fs), dim(
give(_dim)), id(
give(_id)) { }
77 return static_cast<size_t>(fs->size()*dim)==mapper.
mapSize();
108 template <
class E>
struct is_arithmetic{
enum{value=0};};
109 template <>
struct is_arithmetic<real_t>{
enum{value=1};};
110 template <typename E, bool = is_arithmetic<E>::value >
111 class _expr {
using E::GISMO_ERROR_expr;};
113 template<
class T>
class gsFeSpace;
115 template<
class T>
class gsFeSolution;
116 template<
class E>
class symm_expr;
117 template<
class E>
class symmetrize_expr;
118 template<
class E>
class normalized_expr;
119 template<
class E>
class trace_expr;
121 template<
class E>
class adjugate_expr;
122 template<
class E>
class norm_expr;
123 template<
class E>
class sqNorm_expr;
124 template<
class E>
class det_expr;
125 template<
class E>
class value_expr;
127 template<
class E>
class max_expr;
128 template<
class E>
class rowsum_expr;
129 template<
class E>
class colsum_expr;
130 template<
class E>
class col_expr;
131 template<
class T>
class meas_expr;
132 template<
class E>
class inv_expr;
133 template<
class E,
bool cw = false>
class tr_expr;
134 template<
class E>
class cb_expr;
135 template<
class E>
class abs_expr;
141 template<
class T>
class cdiam_expr;
142 template<
class E>
class temp_expr;
143 template<
class E1,
class E2,
bool = E1::ColBlocks && !E1::ScalarValued && !E2::ScalarValued>
class mult_expr
144 {
using E1::GISMO_ERROR_mult_expr_has_invalid_template_arguments;};
148 pow(_expr<E>
const& u, real_t q) {
return pow_expr<E>(u,q); }
153 template <
typename E>
struct expr_traits
157 typedef real_t Scalar;
158 typedef const E Nested_t;
161 # define Temporary_t typename util::conditional<ScalarValued,Scalar, \
162 typename gsMatrix<Scalar>::Base >::type
163 #if __cplusplus >= 201402L || _MSVC_LANG >= 201402L // c++14
164 # define MatExprType auto
165 # define AutoReturn_t auto
167 #else // 199711L, 201103L
168 # define MatExprType typename gsMatrix<real_t>::constRef
169 # define AutoReturn_t typename util::conditional<ScalarValued,real_t,MatExprType>::type
175 template <
typename E>
176 class _expr<E, false>
180 _expr(
const _expr&) { }
190 typedef typename expr_traits<E>::Nested_t Nested_t;
191 typedef typename expr_traits<E>::Scalar Scalar;
197 static_cast<E const&
>(*this).print(os);
216 std::ostream & printDetail(std::ostream &os)
const
218 os << (isVectorTr() ?
"VectorTr " :
219 (isVector() ?
"Vector " :
220 (isMatrix() ?
"Matrix " :
222 <<
"expression of size "<< rows()
223 <<
" x "<<cols()<<
"\n";
230 {
return static_cast<E const&
>(*this).eval(k); }
233 tr_expr<E>
tr()
const
234 {
return tr_expr<E,false>(
static_cast<E const&
>(*this)); }
238 {
return tr_expr<E,true>(
static_cast<E const&
>(*this)); }
241 cb_expr<E>
cb()
const
242 {
return cb_expr<E>(
static_cast<E const&
>(*this)); }
246 {
return sign_expr<E>(
static_cast<E const&
>(*this), tolerance); }
250 {
return exp_expr<E>(
static_cast<E const&
>(*this)); }
259 mult_expr<real_t, ppart_expr<mult_expr<double,E,false>> ,
false>
260 npart()
const {
return -1* ( -(*this) ).ppart() ; }
264 {
return temp_expr<E>(
static_cast<E const&
>(*this)); }
267 inv_expr<E>
const inv()
const
268 {
return inv_expr<E>(
static_cast<E const&
>(*this)); }
272 {
return trace_expr<E>(
static_cast<E const&
>(*this)); }
275 adjugate_expr<E>
adj()
const
276 {
return adjugate_expr<E>(
static_cast<E const&
>(*this)); }
280 {
return norm_expr<E>(
static_cast<E const&
>(*this)); }
284 {
return normalized_expr<E>(
static_cast<E const&
>(*this)); }
288 {
return det_expr<E>(
static_cast<E const&
>(*this)); }
292 {
return sqNorm_expr<E>(
static_cast<E const&
>(*this)); }
295 mult_expr<E,E,0> (sqr)()
const {
return (*
this)*(*this); }
297 symm_expr<E> symm()
const
298 {
return symm_expr<E>(
static_cast<E const&
>(*this)); }
300 symmetrize_expr<E> symmetrize()
const
301 {
return symmetrize_expr<E>(
static_cast<E const&
>(*this)); }
306 {
return value_expr<E>(
static_cast<E const&
>(*this)); }
314 {
return max_expr<E>(
static_cast<E const&
>(*this)); }
318 {
return rowsum_expr<E>(
static_cast<E const&
>(*this)); }
322 {
return colsum_expr<E>(
static_cast<E const&
>(*this)); }
324 col_expr<E> operator[](
const index_t i)
const
325 {
return col_expr<E>(
static_cast<E const&
>(*this),i); }
329 {
return static_cast<E const&
>(*this).rows(); }
333 {
return static_cast<E const&
>(*this).cols(); }
336 {
return static_cast<E const&
>(*this).cardinality_impl(); }
338 static index_t cardinality_impl() {
return 1; }
345 static bool isVectorTr() {
return 2==E::Space; }
346 static bool isMatrix () {
return 3==E::Space; }
351 {
static_cast<E const&
>(*this).parse(evList); }
353 template<
class op>
void apply(op & _op)
const
354 {
static_cast<E const&
>(*this).apply(_op); }
361 return static_cast<E const&
>(*this).rowVar();
369 return static_cast<E const&
>(*this).colVar();
374 operator E&() {
return static_cast< E&
>(*this); }
375 operator E
const&()
const {
return static_cast<const E&
>(*this); }
377 E
const & derived()
const {
return static_cast<const E&
>(*this); }
381 template <
typename E>
382 std::ostream &operator<<(std::ostream &os, const _expr<E> & b)
383 {b.print(os);
return os; }
399 class gsNullExpr :
public _expr<gsNullExpr<T> >
403 operator const gsFeSpace<T> & ()
const
405 static gsFeSpace<T> vv(-1);
413 void parse(gsExprHelper<T> &)
const { }
415 const gsFeSpace<T> & rowVar()
const {
GISMO_ERROR(
"gsNullExpr"); }
416 const gsFeSpace<T> & colVar()
const {
GISMO_ERROR(
"gsNullExpr"); }
418 void print(std::ostream &os)
const { os <<
"NullExpr"; }
420 static const gsNullExpr &
get()
431 class symbol_expr :
public _expr<E>
434 typedef typename expr_traits<E>::Scalar Scalar;
438 const gsFunctionSet<Scalar> * m_fs;
439 const gsFuncData<Scalar> * m_fd;
446 bool isAcross()
const {
return m_isAcross; }
450 E ac(this->derived());
452 ac.m_isAcross =
true;
458 E ac(this->derived());
460 ac.m_isAcross =
false;
465 const gsFunctionSet<Scalar> & source()
const {
return *m_fs;}
468 const gsFuncData<Scalar> & data()
const
470 GISMO_ASSERT(NULL!=m_fd,
"FuncData member not registered "<<
this<<
"/"<< m_fs);
475 void parse(gsExprHelper<Scalar> & evList)
const
481 index_t cardinality_impl()
const
483 GISMO_ASSERT(this->data().actives.rows()!=0,
"Cardinality depends on the NEED_ACTIVE flag");
484 return m_d * this->data().actives.rows();
488 void setSource(
const gsFunctionSet<Scalar> & fs) { m_fs = &fs;}
491 void setData(
const gsFuncData<Scalar> & val) { m_fd = &val;}
492 void setDim(
index_t _d) { m_d = _d; }
493 void clear() { m_fs = NULL; }
496 explicit symbol_expr(
index_t _d)
497 : m_fs(NULL), m_fd(NULL), m_d(_d), m_isAcross(false) { }
500 bool isValid()
const {
return NULL!=m_fd && NULL!=m_fs; }
509 {
return m_fd->values[0].col(k).blockDiag(m_d); }
512 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
513 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
517 GISMO_ASSERT(NULL!=m_fs,
"FeVariable: Function source not registered");
518 return m_fs->targetDim();
521 index_t cols()
const {
return m_d; }
523 void print(std::ostream &os)
const { os <<
"u"; }
528 index_t dim()
const {
return m_d;}
532 index_t targetDim()
const {
return m_fs->targetDim(); }
535 index_t parDim()
const {
return m_fs->domainDim(); }
539 GISMO_ASSERT(0!=m_fd->values[0].size(),
"Probable error.");
540 return m_fd->values[0].rows();
548 class col_expr :
public _expr<col_expr<E> >
550 typename E::Nested_t _c;
553 typedef typename E::Scalar Scalar;
554 typedef const col_expr<E> Nested_t;
556 enum { Space = E::Space, ScalarValued = 0, ColBlocks = 0 };
558 col_expr(
const E & c,
const index_t i) : _c(c), _i(i) { }
565 index_t rows()
const {
return _c.rows(); }
566 index_t cols()
const {
return 1; }
567 void parse(gsExprHelper<Scalar> & evList)
const { _c.parse(evList); }
569 const gsFeSpace<Scalar> & rowVar()
const {
return _c.rowVar(); }
570 const gsFeSpace<Scalar> & colVar()
const {
return _c.colVar(); }
572 void print(std::ostream &os)
const { os<<_c<<
"["<<_i<<
"]"; }
579 class _expr<T, true> :
public _expr<_expr<T> >
584 typedef const _expr<T> Nested_t;
586 explicit _expr(Scalar c) : _c(
give(c)) { }
589 enum {Space = 0, ScalarValued = 1, ColBlocks= 0};
591 inline Scalar eval(
const index_t )
const {
return _c; }
593 inline _expr val()
const {
return *
this; }
594 index_t rows()
const {
return 0; }
595 index_t cols()
const {
return 0; }
596 void parse(gsExprHelper<Scalar> &)
const { }
598 const gsFeSpace<T> & rowVar()
const {
return gsNullExpr<T>::get(); }
599 const gsFeSpace<T> & colVar()
const {
return gsNullExpr<T>::get(); }
601 void print(std::ostream &os)
const { os<<_c; }
608 class gsGeometryMap :
public _expr<gsGeometryMap<T> >
610 const gsFunctionSet<T> * m_fs;
611 const gsMapData<T> * m_fd;
617 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
619 bool isAcross()
const {
return m_isAcross; }
621 gsGeometryMap right()
const
625 ac.m_isAcross =
true;
629 gsGeometryMap left()
const
633 ac.m_isAcross =
false;
638 const gsFunctionSet<T> & source()
const {
return *m_fs;}
641 const gsMapData<T> & data()
const
643 GISMO_ASSERT(NULL!=m_fd,
"gsGeometryMap: invalid data "<< m_fs <<
","<<m_fd);
647 index_t targetDim()
const {
return m_fs->targetDim();}
648 index_t domainDim()
const {
return m_fs->domainDim();}
651 void copyCoefs(
const gsGeometryMap<T> & other)
const
653 GISMO_ASSERT(
dynamic_cast<const gsMultiPatch<T>*
>( this->m_fs ),
"error");
654 const gsMultiPatch<T> & thisMP =
static_cast<const gsMultiPatch<T>&
>(*this->m_fs );
655 GISMO_ASSERT(
dynamic_cast<const gsMultiPatch<T>*
>( other.m_fs ),
"error");
656 const gsMultiPatch<T> & otherMP =
static_cast<const gsMultiPatch<T>&
>(*other.m_fs );
657 GISMO_ASSERT( (thisMP.domainDim()==otherMP.domainDim())&&
658 (thisMP.geoDim()==otherMP.geoDim())&&
659 (thisMP.coefsSize() == otherMP.coefsSize())&&
660 (thisMP.nPatches()==otherMP.nPatches()),
661 "The geometryMaps are not compatible!");
664 for (
index_t p=0; p < thisMP.nPatches(); p++ )
667 thisMP.patch(p).coefs() = otherMP.patch(p).coefs();
672 void deformBy(
const gsFeSolution<T> & displacement)
const
674 const index_t dim = m_fs->domainDim();
676 const gsMultiBasis<T> & mb =
static_cast<const gsMultiBasis<T>&
>(displacement.space().source());
677 const gsMultiPatch<T> & mp =
static_cast<const gsMultiPatch<T>&
>(*this->m_fs );
678 GISMO_ASSERT(
dynamic_cast<const gsMultiBasis<T>*
>(&displacement.space().source()),
"error");
679 GISMO_ASSERT(
dynamic_cast<const gsMultiPatch<T>*
>( this->m_fs),
"error");
682 for (
size_t p=0; p < mp.nPatches(); p++ )
685 gsMatrix<T> &result = mp.patch(p).coefs();
688 const index_t sz = mb[p].size();
691 for (
index_t c = 0; c!=dim; c++)
694 for (
index_t i = 0; i < sz; ++i)
696 const int ii = displacement.mapper().index(i, p, c);
697 if ( displacement.mapper().is_free_index(ii) )
699 result(i,c) += displacement.coefs().at(ii);
703 result(i,c) += displacement.fixedPart().at( displacement.mapper().global_to_bindex(ii));
715 void print(std::ostream &os)
const { os <<
"G"; }
717 auto eval(
const index_t k)
const -> decltype(m_fd->values[0].col(k))
718 {
return m_fd->values[0].col(k); }
722 gsGeometryMap() : m_fs(NULL), m_fd(NULL), m_isAcross(false) { }
724 void setSource(
const gsFunctionSet<Scalar> & fs) { m_fs = &fs;}
725 void setData(
const gsMapData<Scalar> & val) { m_fd = &val;}
729 index_t rows()
const {
return m_fs->targetDim(); }
730 index_t cols()
const {
return 1; }
732 const gsFeSpace<T> & rowVar()
const {
return gsNullExpr<T>::get(); }
733 const gsFeSpace<T> & colVar()
const {
return gsNullExpr<T>::get(); }
735 void parse(gsExprHelper<Scalar> & evList)
const
743 template <
typename T>
struct expr_traits<gsGeometryMap<T> >
746 typedef const gsGeometryMap<T> Nested_t;
753 class meas_expr :
public _expr<meas_expr<T> >
755 typename gsGeometryMap<T>::Nested_t _G;
758 enum {Space = 0, ScalarValued = 1, ColBlocks = 0};
762 meas_expr(
const gsGeometryMap<T> & G) : _G(G) { }
766 return _G.data().measures.at(k);
769 index_t rows()
const {
return 0; }
770 index_t cols()
const {
return 0; }
771 void parse(gsExprHelper<Scalar> & evList)
const
777 const gsFeSpace<T> & rowVar()
const {
return gsNullExpr<T>::get(); }
778 const gsFeSpace<T> & colVar()
const {
return gsNullExpr<T>::get(); }
780 void print(std::ostream &os)
const { os <<
"meas("; _G.print(os); os <<
")"; }
790 friend class cdiam_expr<T>;
792 const gsDomainIterator<T> * m_di;
794 const gsVector<T> * m_weights;
797 gsFeElement(
const gsFeElement &);
801 gsFeElement() : m_di(NULL), m_weights(nullptr) { }
803 void set(
const gsDomainIterator<T> & di,
const gsVector<T> & weights)
804 { m_di = &di, m_weights = &weights; }
806 bool isValid()
const {
return nullptr!=m_weights; }
808 const gsVector<T> & weights()
const {
return *m_weights;}
811 integral_expr<E> integral(
const _expr<E>& ff)
const
812 {
return integral_expr<E>(*
this,ff); }
814 typedef integral_expr<T> AreaRetType;
815 AreaRetType area()
const
816 {
return integral(_expr<T,true>(1)); }
818 typedef integral_expr<meas_expr<T> > PHAreaRetType;
820 PHAreaRetType area(
const gsGeometryMap<Scalar> & _G)
const
821 {
return integral(meas_expr<T>(_G)); }
823 typedef pow_expr<integral_expr<T> > DiamRetType;
825 DiamRetType diam() const
826 {
return pow(integral(_expr<T,true>(1)),(T)(1)/(T)(2)); }
828 typedef pow_expr<integral_expr<meas_expr<T> > > PHDiamRetType;
830 PHDiamRetType diam(
const gsGeometryMap<Scalar> & _G)
const
831 {
return pow(integral(meas_expr<T>(_G)),(T)(1)/(T)(2)); }
837 void print(std::ostream &os)
const { os <<
"e"; }
839 void parse(gsExprHelper<T> & evList)
const
849 class integral_expr :
public _expr<integral_expr<E> >
853 typedef real_t Scalar;
854 mutable Scalar m_val;
856 const gsFeElement<Scalar> &
_e;
857 typename _expr<E>::Nested_t _ff;
859 enum {Space= 0, ScalarValued= 1, ColBlocks = 0};
861 integral_expr(
const gsFeElement<Scalar> & el,
const _expr<E> & u)
862 : m_val(-1),
_e(el), _ff(u) { }
864 const Scalar & eval(
const index_t k)
const
866 GISMO_ENSURE(
_e.isValid(),
"Element is valid within integrals only.");
869 const Scalar * w =
_e.weights().data();
870 m_val = (*w) * _ff.val().eval(0);
871 for (
index_t j = 1; j !=
_e.weights().rows(); ++j)
872 m_val += (*(++w)) * _ff.val().eval(j);
877 inline const integral_expr<E> & val()
const {
return *
this; }
878 inline index_t rows()
const {
return 0; }
879 inline index_t cols()
const {
return 0; }
880 void parse(gsExprHelper<Scalar> & evList)
const
885 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get(); }
886 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get(); }
888 void print(std::ostream &os)
const
922 template <
typename T>
923 struct expr_traits<gsFeVariable<T> >
926 typedef const gsFeVariable<T> Nested_t;
934 class gsFeVariable :
public symbol_expr< gsFeVariable<T> >
937 typedef symbol_expr< gsFeVariable<T> > Base;
939 explicit gsFeVariable(
index_t _d = 1) : Base(_d) { }
941 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
945 class gsComposition :
public symbol_expr< gsComposition<T> >
948 typedef symbol_expr< gsComposition<T> > Base;
949 typename gsGeometryMap<T>::Nested_t _G;
951 explicit gsComposition(
const gsGeometryMap<T> & G,
index_t _d = 1)
952 : Base(_d), _G(G) { }
954 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
956 typename gsMatrix<T>::constColumn
957 eval(
const index_t k)
const {
return this->m_fd->values[0].col(k); }
959 const gsGeometryMap<T> & inner()
const {
return _G;};
961 void parse(gsExprHelper<T> & evList)
const
979 class gsFeSpace :
public symbol_expr< gsFeSpace<T> >
981 friend class gsNullExpr<T>;
984 typedef symbol_expr< gsFeSpace<T> > Base;
987 gsFeSpaceData<T> * m_sd;
990 enum{Space = 1, ScalarValued=0, ColBlocks=0};
992 typedef const gsFeSpace Nested_t;
996 const gsFeSpace<T> & rowVar()
const {
return *
this;}
998 gsDofMapper & mapper()
1000 GISMO_ASSERT(NULL!=m_sd,
"Space/mapper not properly initialized.");
1001 return m_sd->mapper;
1004 const gsDofMapper & mapper()
const
1005 {
return const_cast<gsFeSpace*
>(
this)->mapper();}
1007 inline const gsMatrix<T> & fixedPart()
const {
return m_sd->fixedDofs;}
1008 gsMatrix<T> & fixedPart() {
return m_sd->fixedDofs;}
1010 index_t id()
const {
return (m_sd ? m_sd->id : -101); }
1011 void setSpaceData(gsFeSpaceData<T>& sd) {m_sd = &sd;}
1013 index_t interfaceCont()
const {
return m_sd->cont;}
1016 GISMO_ASSERT(_r>-2 && _r<1,
"Invalid or not implemented (r="<<_r<<
").");
1017 return m_sd->cont = _r;
1020 gsFeSolution<T>
function(
const gsMatrix<T>& solVector)
const
1021 {
return gsFeSolution<T>(*this); }
1023 void getCoeffs(
const gsMatrix<T>& solVector, gsMatrix<T> & result,
1026 const index_t dim = this->dim();
1028 const gsMultiBasis<T> & mb =
static_cast<const gsMultiBasis<T>&
>(this->source());
1029 GISMO_ASSERT(
dynamic_cast<const gsMultiBasis<T>*
>(&this->source()),
"error");
1032 const index_t sz = mb[p].size();
1033 result.resize(sz, dim!=1 ? dim : solVector.cols());
1035 for (
index_t c = 0; c!=dim; c++)
1038 for (
index_t i = 0; i < sz; ++i)
1040 const int ii = m_sd->mapper.index(i, p, c);
1041 if ( m_sd->mapper.is_free_index(ii) )
1042 for(
index_t s = 0; s != solVector.cols(); ++s )
1043 result(i,c+s) = solVector(ii,s);
1045 result(i,c) = m_sd->fixedDofs.at( m_sd->mapper.global_to_bindex(ii) );
1053 void setupMapper(gsDofMapper dofsMapper)
const
1055 GISMO_ASSERT( dofsMapper.isFinalized(),
"The provided dof-mapper is not finalized.");
1056 GISMO_ASSERT( dofsMapper.mapSize()==
static_cast<size_t>(this->source().size()*dofsMapper.numComponents()),
"The dof-mapper is not consistent: mapSize()="<<dofsMapper.mapSize()<<
"!="<<
static_cast<size_t>(this->source().size())<<
"=this->source().size()");
1057 m_sd->mapper =
give(dofsMapper);
1060 void setup(
const index_t _icont = -1)
const
1062 this->setInterfaceCont(_icont);
1063 m_sd->mapper = gsDofMapper();
1065 if (
const gsMultiBasis<T> * mb =
1066 dynamic_cast<const gsMultiBasis<T>*
>(&this->source()) )
1068 m_sd->mapper = gsDofMapper(*mb, this->dim() );
1070 if ( 0==this->interfaceCont() )
1072 for ( gsBoxTopology::const_iiterator it = mb->topology().iBegin();
1073 it != mb->topology().iEnd(); ++it )
1075 mb->matchInterface(*it, m_sd->mapper);
1080 if (
const gsMappedBasis<2,T> * mb =
1081 dynamic_cast<const gsMappedBasis<2,T>*
>(&this->source()) )
1083 m_sd->mapper.setIdentity(mb->nPatches(), mb->size() , this->dim());
1086 m_sd->mapper.finalize();
1089 void setup(
const gsBoundaryConditions<T> & bc,
const index_t dir_values,
1090 const index_t _icont = -1)
const
1092 this->setInterfaceCont(_icont);
1093 m_sd->mapper = gsDofMapper();
1094 const index_t dim = this->dim();
1095 const gsMultiBasis<T> *mb =
dynamic_cast<const gsMultiBasis<T> *
>(&this->source());
1098 m_sd->mapper = gsDofMapper(*mb, this->dim());
1100 if (0 == this->interfaceCont())
1102 for (gsBoxTopology::const_iiterator it = mb->topology().iBegin();
1103 it != mb->topology().iEnd(); ++it) {
1104 mb->matchInterface(*it, m_sd->mapper);
1109 gsMatrix<index_t> bnd;
1110 for (
typename gsBoundaryConditions<T>::const_iterator
1111 it = bc.begin(
"Dirichlet"); it != bc.end(
"Dirichlet"); ++it)
1113 const index_t cc = it->unkComponent();
1114 GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1115 "Problem: a boundary condition is set on a patch id which does not exist.");
1117 bnd = mb->basis(it->ps.patch).boundary(it->ps.side());
1118 m_sd->mapper.markBoundary(it->ps.patch, bnd, cc);
1121 gsMatrix<index_t> bnd1;
1122 for (
typename gsBoundaryConditions<T>::const_iterator
1123 it = bc.begin(
"Clamped"); it != bc.end(
"Clamped"); ++it)
1125 const index_t cc = it->unkComponent();
1127 GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1128 "Problem: a boundary condition is set on a patch id which does not exist.");
1130 bnd = mb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 0);
1131 bnd1 = mb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 1);
1133 if (!it->ps.parameter())
1135 for (
index_t c = 0; c!=dim; c++)
1137 if (c==cc || cc==-1 )
1138 for (
index_t k = 0; k < bnd.size(); ++k)
1139 m_sd->mapper.matchDof(it->ps.patch, (bnd)(k, 0),
1140 it->ps.patch, (bnd1)(k, 0), c);
1146 for (
typename gsBoundaryConditions<T>::const_iterator
1147 it = bc.begin(
"Collapsed"); it != bc.end(
"Collapsed"); ++it)
1149 const index_t cc = it->unkComponent();
1151 GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1152 "Problem: a boundary condition is set on a patch id which does not exist.");
1154 bnd = mb->basis(it->ps.patch).boundary(it->ps.side());
1157 for (
index_t c = 0; c!=dim; c++)
1159 if (c==cc || cc==-1)
1160 for (
index_t k = 0; k < bnd.size() - 1; ++k)
1161 m_sd->mapper.matchDof(it->ps.patch, (bnd)(0, 0),
1162 it->ps.patch, (bnd)(k + 1, 0), c);
1167 for (
typename gsBoundaryConditions<T>::const_cpliterator
1168 it = bc.coupledBegin(); it != bc.coupledEnd(); ++it)
1170 const index_t cc = it->component;
1173 "Problem: a boundary condition is set on a patch id which does not exist.");
1175 "Problem: a boundary condition is set on a patch id which does not exist.");
1178 bnd = mb->basis(it->ifc.first().patch).boundary(it->ifc.first().side());
1179 bnd1 = mb->basis(it->ifc.second().patch).boundary(it->ifc.second().side());
1182 for (
index_t c = 0; c!=dim; c++)
1184 if (c==cc || cc==-1)
1186 for (
index_t k = 0; k < bnd.size() - 1; ++k)
1187 m_sd->mapper.matchDof(it->ifc.first() .patch, (bnd)(0, 0),
1188 it->ifc.first() .patch, (bnd)(k + 1, 0), c);
1189 for (
index_t k = 0; k < bnd1.size(); ++k)
1190 m_sd->mapper.matchDof(it->ifc.first() .patch, (bnd)(0, 0),
1191 it->ifc.second().patch, (bnd1)(k, 0), c);
1197 for (
typename gsBoundaryConditions<T>::const_citerator
1198 it = bc.cornerBegin(); it != bc.cornerEnd(); ++it)
1200 for (
index_t r = 0; r!=this->dim(); ++r)
1202 if (it->component!=-1 && r!=it->component)
continue;
1205 GISMO_ASSERT(static_cast<size_t>(it->patch) < mb->nBases(),
1206 "Problem: a corner boundary condition is set on a patch id which does not exist.");
1207 m_sd->mapper.eliminateDof(mb->basis(it->patch).functionAtCorner(it->corner),
1208 it->patch, it->component);
1212 }
else if (
const gsBasis<T> *b =
1213 dynamic_cast<const gsBasis<T> *
>(&this->source()))
1215 m_sd->mapper = gsDofMapper(*b, this->dim() );
1216 gsMatrix<index_t> bnd;
1217 for (
typename gsBoundaryConditions<T>::const_iterator
1218 it = bc.begin(
"Dirichlet"); it != bc.end(
"Dirichlet"); ++it) {
1220 "Problem: a boundary condition is set on a patch id which does not exist.");
1222 bnd = b->boundary(it->ps.side());
1223 m_sd->mapper.markBoundary(0, bnd, it->unkComponent());
1226 for (
typename gsBoundaryConditions<T>::const_iterator
1227 it = bc.begin(
"Clamped"); it != bc.end(
"Clamped"); ++it) {
1229 "Problem: a boundary condition is set on a patch id which does not exist.");
1231 bnd = b->boundary(it->ps.side());
1236 m_sd->mapper = gsDofMapper(*b);
1237 for (
typename gsBoundaryConditions<T>::const_iterator
1238 it = bc.begin(
"Collapsed"); it != bc.end(
"Collapsed"); ++it) {
1240 "Problem: a boundary condition is set on a patch id which does not exist.");
1242 bnd = b->boundary(it->ps.side());
1246 }
else if (
const gsMappedBasis<2, T> *mapb =
1247 dynamic_cast<const gsMappedBasis<2, T> *
>(&this->source()))
1249 m_sd->mapper.setIdentity(mapb->nPatches(), mapb->size(), this->dim());
1251 if (0 == this->interfaceCont())
1253 gsMatrix<index_t> int1, int2;
1254 for (gsBoxTopology::const_iiterator it = mapb->getTopol().iBegin();
1255 it != mapb->getTopol().iEnd(); ++it) {
1256 int1 = mapb->basis(it->first().patch).boundaryOffset(it->first().side(), 0);
1257 int2 = mapb->basis(it->second().patch).boundaryOffset(it->second().side(), 0);
1259 m_sd->mapper.matchDofs(it->first().patch, int1, it->second().patch, int2);
1262 if (1 == this->interfaceCont())
1264 GISMO_ERROR(
"Boundary offset function is not implemented for gsMappedBasis in general.");
1267 gsMatrix<index_t> bnd;
1268 for (
typename gsBoundaryConditions<T>::const_iterator
1269 it = bc.begin(
"Dirichlet"); it != bc.end(
"Dirichlet"); ++it) {
1270 const index_t cc = it->unkComponent();
1271 GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1272 "Problem: a boundary condition is set on a patch id which does not exist.");
1274 bnd = mapb->basis(it->ps.patch).boundary(it->ps.side());
1275 m_sd->mapper.markBoundary(it->ps.patch, bnd, cc);
1279 gsMatrix<index_t> bnd1;
1280 for (
typename gsBoundaryConditions<T>::const_iterator
1281 it = bc.begin(
"Clamped"); it != bc.end(
"Clamped"); ++it) {
1282 const index_t cc = it->unkComponent();
1284 GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1285 "Problem: a boundary condition is set on a patch id which does not exist.");
1287 bnd = mapb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 0);
1288 bnd1 = mapb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 1);
1293 if (!it->ps.parameter())
1295 for (
index_t c = 0; c!=dim; c++)
1297 if (c==cc || cc==-1 )
1298 for (
index_t k = 0; k < bnd.size() - 1; ++k)
1299 m_sd->mapper.matchDof( it->ps.patch, (bnd)(k, 0),
1300 it->ps.patch, (bnd1)(k, 0), c);
1303 gsWarn <<
"Unable to apply clamped condition.\n";
1307 for (
typename gsBoundaryConditions<T>::const_iterator
1308 it = bc.begin(
"Collapsed"); it != bc.end(
"Collapsed"); ++it) {
1309 const index_t cc = it->unkComponent();
1311 GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1312 "Problem: a boundary condition is set on a patch id which does not exist.");
1314 bnd = mapb->basis(it->ps.patch).boundary(it->ps.side());
1320 for (
index_t c = 0; c!=dim; c++)
1322 if (c==cc || cc==-1)
1323 for (
index_t k = 0; k < bnd.size() - 1; ++k)
1324 m_sd->mapper.matchDof(it->ps.patch, (bnd)(0, 0),
1325 it->ps.patch, (bnd)(k + 1, 0), c);
1331 for (
typename gsBoundaryConditions<T>::const_citerator
1332 it = bc.cornerBegin(); it != bc.cornerEnd(); ++it)
1336 "Problem: a corner boundary condition is set on a patch id which does not exist.");
1337 m_sd->mapper.eliminateDof(mapb->basis(it->patch).functionAtCorner(it->corner), it->patch, it->component);
1341 GISMO_ASSERT(0 == bc.size(),
"Problem: BCs are ignored.");
1342 m_sd->mapper.setIdentity(this->source().nPieces(), this->source().size());
1345 m_sd->mapper.finalize();
1348 gsDirichletValues(bc, dir_values, *
this);
1351 void print(std::ostream &os)
const { os <<
"u"; }
1355 friend class symbol_expr<gsFeSpace>;
1356 explicit gsFeSpace(
index_t _d = 1) : Base(_d), m_sd(nullptr) { }
1359 template<
class T>
inline bool
1360 operator== (
const gsFeSpace<T> & a,
const gsFeSpace<T> & b)
1361 {
return a.id()== b.id() && a.isAcross()==b.isAcross(); }
1371 class gsFeSolution :
public _expr<gsFeSolution<T> >
1374 const gsFeSpace<T> _u;
1380 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
1382 bool isAcross()
const {
return m_isAcross; }
1384 gsFeSolution right()
const
1386 gsFeSolution ac(*
this);
1387 ac.m_isAcross =
true;
1391 gsFeSolution left()
const {
return gsFeSolution(*
this); }
1393 explicit gsFeSolution(
const gsFeSpace<T> & u) : _u(u), _Sv(NULL) { }
1395 gsFeSolution(
const gsFeSpace<T> & u, gsMatrix<T> & Sv) : _u(u), _Sv(&Sv) { }
1397 const gsFeSpace<T> & space()
const {
return _u;};
1399 mutable gsMatrix<T> res;
1400 const gsMatrix<T> & eval(
index_t k)
const
1402 bool singleActives = (1 == _u.data().actives.cols());
1404 res.setZero(_u.dim(), 1);
1405 const gsDofMapper & map = _u.mapper();
1406 GISMO_ASSERT(_Sv->size()==map.freeSize(),
"The solution vector has wrong dimensions: "<<_Sv->size()<<
" != "<<map.freeSize());
1408 for (
index_t c = 0; c!=_u.dim(); c++)
1410 for (
index_t i = 0; i!=_u.data().actives.rows(); ++i)
1412 const index_t ii = map.index(_u.data().actives(i, singleActives ? 0 : k), _u.data().patchId, c);
1413 if ( map.is_free_index(ii) )
1414 res.at(c) += _Sv->at(ii) * _u.data().values[0](i,k);
1416 res.at(c) += _u.data().values[0](i,k) *
1417 _u.fixedPart().at( map.global_to_bindex(ii) );
1427 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
1428 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
1429 index_t rows()
const {
return _u.dim(); }
1431 static index_t cols() {
return 1; }
1433 void parse(gsExprHelper<Scalar> & evList)
const
1439 void print(std::ostream &os)
const { os <<
"s"; }
1442 index_t dim()
const {
return _u.dim();}
1445 {
return _u.source().domainDim(); }
1448 const gsDofMapper & mapper()
const {
return _u.mapper();}
1450 inline const gsMatrix<T> & fixedPart()
const {
return _u.fixedPart();}
1451 gsMatrix<T> & fixedPart() {
return _u.fixedPart();}
1454 const gsFuncData<T> & data()
const {
return _u.data();}
1456 void setSolutionVector(gsMatrix<T>& solVector)
1457 { _Sv = & solVector; }
1464 void setComponent(
index_t component, real_t value,
index_t patch=-1)
1466 gsMatrix<T> & solVector = *_Sv;
1467 const gsDofMapper & mapper = _u.mapper();
1472 patchEnd = _u.mapper().numPatches();
1476 patchEnd = patch + 1;
1479 for (
index_t p=patchStart; p!=patchEnd; ++p)
1481 for (
size_t i = 0; i != mapper.patchSize(p, component); ++i)
1483 const index_t ii = mapper.index(i, p, component);
1484 if ( mapper.is_free_index(ii) )
1485 solVector.at(ii) = value;
1490 const gsMatrix<T> & coefs()
const {
return *_Sv; }
1501 GISMO_ASSERT(j<_Sv->size(),
"Solution vector is not initialized/allocated, sz="<<_Sv->size() );
1509 void extract(gsMatrix<T> & result,
const index_t p = 0)
const
1510 { _u.getCoeffs(*_Sv, result, p); }
1514 void extractFull(gsMatrix<T> & result)
const
1518 const size_t totalSz = _u.mapper().mapSize();
1519 result.resize(totalSz, 1);
1520 for (
size_t p=0; p!=_u.mapper().numPatches(); ++p)
1522 offset = _u.mapper().offset(p);
1525 for (
index_t c = 0; c!=dim; c++)
1527 const index_t sz = _u.mapper().patchSize(p,c);
1530 for (
index_t i = 0; i < sz; ++i)
1533 const int ii = _u.mapper().index(i, p, c);
1535 if ( _u.mapper().is_free_index(ii) )
1537 result(i+offset,0) = _Sv->at(ii);
1540 result(i+offset,0) = _u.fixedPart().at( _u.mapper().global_to_bindex(ii) );
1548 void extract(gsMultiPatch<T> & result)
const
1552 if(
const gsMultiBasis<T>* basis =
dynamic_cast<const gsMultiBasis<T>*
>(&_u.source()) )
1553 for (
size_t i = 0; i != basis->nBases(); ++i)
1555 memory::unique_ptr<gsGeometry<T> > p(this->extractPiece(i));
1556 result.addPatch(*p);
1561 void extract(gsMappedSpline<2,T> & result)
const
1563 if(
const gsMappedBasis<2,T>* basis =
dynamic_cast<const gsMappedBasis<2,T>*
>(&_u.source()) )
1566 this->extractFull(coefs);
1567 coefs.resize(coefs.rows()/_u.dim(),_u.dim());
1568 result.init(*basis,coefs);
1573 memory::unique_ptr<gsGeometry<T> > extractPiece(
const index_t p)
const
1575 if (
const gsBasis<T> * b =
dynamic_cast<const gsBasis<T>*
>(&_u.source().piece(p)) )
1579 return b->makeGeometry(
give(cf));
1585 void insert(
const gsGeometry<T> & g,
const index_t p = 0)
const
1587 const gsMatrix<T> & cf = g.coefs();
1588 gsMatrix<T> & sol = *_Sv;
1590 const gsDofMapper & mapper = _u.mapper();
1591 for (
index_t c = 0; c!=_u.dim(); c++)
1593 for (
index_t i = 0; i != cf.rows(); ++i)
1595 const index_t ii = mapper.index(i, p, c);
1596 if ( mapper.is_free_index(ii) )
1597 sol.at(ii) = cf(i, c);
1612 template<
class E,
bool cw>
1613 class tr_expr :
public _expr<tr_expr<E,cw> >
1615 typename E::Nested_t _u;
1619 typedef typename E::Scalar Scalar;
1621 tr_expr(_expr<E>
const& u)
1625 enum {ColBlocks = E::ColBlocks, ScalarValued=E::ScalarValued};
1626 enum {Space = cw?E::Space:(E::Space==1?2:(E::Space==2?1:E::Space))};
1628 mutable Temporary_t res;
1629 const Temporary_t & eval(
const index_t k)
const
1632 res = _u.eval(k).blockTranspose( _u.cardinality() );
1634 res = _u.eval(k).transpose();
1638 index_t rows()
const {
return _u.cols(); }
1640 index_t cols()
const {
return _u.rows(); }
1642 void parse(gsExprHelper<Scalar> & evList)
const
1643 { _u.parse(evList); }
1645 const gsFeSpace<Scalar> & rowVar()
const {
return cw?_u.rowVar():_u.colVar(); }
1646 const gsFeSpace<Scalar> & colVar()
const {
return cw?_u.colVar():_u.rowVar(); }
1648 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
1650 void print(std::ostream &os)
const { os<<
"("; _u.print(os); os <<
")\u1D40"; }
1667 class cb_expr :
public _expr<cb_expr<E> >
1669 typename E::Nested_t _u;
1673 typedef typename E::Scalar Scalar;
1675 cb_expr(_expr<E>
const& u)
1679 enum {ColBlocks = 1, ScalarValued=E::ScalarValued};
1680 enum {Space = E::Space};
1682 mutable gsMatrix<Scalar> ev, res;
1684 const gsMatrix<Scalar> & eval(
const index_t k)
const
1697 GISMO_ASSERT(res.rows() % _u.rows() == 0 && res.cols() % _u.cols() == 0,
"Result is not a multiple of the space dimensions?");
1700 if ( (cardinality = res.rows() / _u.rows()) >= 1 && res.cols() / _u.cols() == 1 )
1702 res.resize(_u.rows(), cardinality * _u.cols());
1703 for (
index_t r = 0; r!=cardinality; r++)
1704 res.block(0 , r * _u.cols(), _u.rows(), _u.cols()) = ev.block( r * _u.rows(), 0, _u.rows(), _u.cols() );
1706 else if ( (cardinality = res.rows() / _u.rows()) == 1 && res.cols() / _u.cols() >= 1 )
1708 res.resize(_u.rows(), cardinality * _u.cols());
1709 for (
index_t r = 0; r!=cardinality; r++)
1710 res.block(0 , r * _u.cols(), _u.rows(), _u.cols()) = ev.block( 0, r * _u.cols(), _u.rows(), _u.cols() );
1716 index_t rows()
const {
return _u.rows(); }
1718 index_t cols()
const {
return _u.cols(); }
1720 void parse(gsExprHelper<Scalar> & evList)
const
1721 { _u.parse(evList); }
1723 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
1724 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
1726 index_t cardinality_impl()
const
1730 if ( res.rows() / _u.rows() >= 1 && res.cols() / _u.cols() == 1 )
1731 cardinality = res.rows() / _u.rows();
1732 else if ( res.rows() / _u.rows() == 1 && res.cols() / _u.cols() >= 1 )
1733 cardinality = res.cols() / _u.cols();
1735 GISMO_ERROR(
"Cardinality for cb_expr cannot be determined.");
1740 void print(std::ostream &os)
const { os<<
"{"; _u.print(os); os <<
"}"; }
1748 class temp_expr :
public _expr<temp_expr<E> >
1750 typename E::Nested_t _u;
1751 typedef typename E::Scalar Scalar;
1752 mutable gsMatrix<Scalar> tmp;
1755 temp_expr(_expr<E>
const& u)
1759 enum {Space = E::Space, ScalarValued = E::ScalarValued,
1760 ColBlocks = E::ColBlocks};
1764 const gsMatrix<Scalar> & eval(
const index_t k)
const
1770 index_t rows()
const {
return _u.rows(); }
1771 index_t cols()
const {
return _u.cols(); }
1772 void parse(gsExprHelper<Scalar> & evList)
const { _u.parse(evList); }
1773 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
1774 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
1776 void print(std::ostream &os)
const { _u.print(os); }
1783 class trace_expr :
public _expr<trace_expr<E> >
1786 typedef typename E::Scalar Scalar;
1787 enum {ScalarValued = 0, Space = E::Space, ColBlocks= 0};
1790 typename E::Nested_t _u;
1791 mutable gsMatrix<Scalar> res;
1794 trace_expr(_expr<E>
const& u) : _u(u)
1801 const gsMatrix<Scalar> & eval(
const index_t k)
const
1803 auto tmp = _u.eval(k);
1805 const index_t r = _u.cardinality();
1811 for (
index_t i = 0; i!=r; ++i)
1812 res.at(i) = tmp.middleCols(i*cb,cb).trace();
1819 index_t rows()
const {
return _u.cols() / _u.rows(); }
1820 index_t cols()
const {
return 1; }
1822 index_t cardinality_impl()
const {
return _u.cardinality(); }
1824 void parse(gsExprHelper<Scalar> & evList)
const
1825 { _u.parse(evList); }
1827 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
1828 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
1830 void print(std::ostream &os)
const { os <<
"trace("; _u.print(os); os<<
")"; }
1837 class adjugate_expr :
public _expr<adjugate_expr<E> >
1840 typedef typename E::Scalar Scalar;
1841 enum {ScalarValued = 0, ColBlocks = E::ColBlocks};
1842 enum {Space = E::Space};
1844 typename E::Nested_t _u;
1845 mutable gsMatrix<Scalar> res;
1848 adjugate_expr(_expr<E>
const& u) : _u(u)
1855 const gsMatrix<Scalar> & eval(
const index_t k)
const
1857 auto tmp = _u.eval(k);
1859 const index_t r = _u.cols() / cb;
1860 res.resize(_u.rows(),_u.cols());
1861 for (
index_t i = 0; i!=r; ++i){
1862 res.middleCols(i*cb,cb) = tmp.middleCols(i*cb,cb).adjugate();
1870 index_t rows()
const {
return _u.rows(); }
1871 index_t cols()
const {
return _u.cols(); }
1873 void parse(gsExprHelper<Scalar> & evList)
const
1874 { _u.parse(evList); }
1876 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
1877 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
1879 void print(std::ostream &os)
const { os <<
"adj("; _u.print(os); os<<
")"; }
1883 class reshape_expr :
public _expr<reshape_expr<E> >
1886 typedef typename E::Scalar Scalar;
1887 enum {ScalarValued = 0, ColBlocks = E::ColBlocks};
1888 enum {Space = E::Space};
1890 typename E::Nested_t _u;
1892 mutable gsMatrix<Scalar> tmp;
1897 reshape_expr(_expr<E>
const& u,
index_t n,
index_t m) : _u(u), _n(n), _m(m)
1902 const gsAsConstMatrix<Scalar> eval(
const index_t k)
const
1905 GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m,
"Wrong dimension "
1906 << _u.rows() <<
" x "<<_u.cols() <<
"!=" << _n <<
" * "<< _m );
1908 return gsAsConstMatrix<Scalar>(tmp.data(),_n,_m);
1911 index_t rows()
const {
return _n; }
1912 index_t cols()
const {
return _m; }
1914 void parse(gsExprHelper<Scalar> & evList)
const
1915 { _u.parse(evList); }
1917 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
1918 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
1920 void print(std::ostream &os)
const { os <<
"reshape("; _u.print(os); os<<
","<<_n<<
","<<_m<<
")"; }
1924 template <
typename E> EIGEN_STRONG_INLINE
1926 {
return reshape_expr<E>(u, n, m); }
1929 class replicate_expr :
public _expr<replicate_expr<E> >
1932 typedef typename E::Scalar Scalar;
1933 enum {ScalarValued = 0, Space = E::Space, ColBlocks= E::ColBlocks};
1935 typename E::Nested_t _u;
1937 mutable gsMatrix<Scalar> tmp;
1942 replicate_expr(_expr<E>
const& u,
index_t n,
index_t m) : _u(u), _n(n), _m(m)
1946 auto eval(
const index_t k)
const -> decltype(tmp.replicate(_n,_m))
1949 return tmp.replicate(_n,_m);
1952 index_t rows()
const {
return _n*_u.rows(); }
1953 index_t cols()
const {
return _m*_u.cols(); }
1955 void parse(gsExprHelper<Scalar> & evList)
const
1956 { _u.parse(evList); }
1958 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
1959 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
1960 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
1962 void print(std::ostream &os)
const { os <<
"replicate("; _u.print(os); os<<
","<<_n<<
","<<_m<<
")"; }
1966 template <
typename E> EIGEN_STRONG_INLINE
1968 {
return replicate_expr<E>(u, n, m); }
1981 typedef typename E::Scalar Scalar;
1982 enum {ScalarValued = 0, Space = E::Space, ColBlocks= 0};
1984 typename E::Nested_t _u;
1997 const index_t numActives = _u.cardinality();
1999 for (
index_t i = 0; i<numActives; ++i)
2001 tmp(0,2*i+1) += tmp(1,2*i);
2002 std::swap(tmp(1,2*i), tmp(1,2*i+1));
2005 tmp.resize(4,numActives);
2006 tmp.conservativeResize(3,numActives);
2009 tmp.transposeInPlace();
2011 tmp.transposeInPlace();
2016 index_t rows()
const {
return 1; }
2017 index_t cols()
const {
return 3; }
2020 { _u.parse(evList); }
2024 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
2026 void print(std::ostream &os)
const { os <<
"flat("; _u.print(os); os<<
")"; }
2030 template <
typename E> EIGEN_STRONG_INLINE
2039 class diag_expr :
public _expr<diag_expr<E> >
2042 typedef typename E::Scalar Scalar;
2043 enum {Space=0, ColBlocks=E::ColBlocks, ScalarValued = 0};
2045 typename E::Nested_t _u;
2046 mutable gsMatrix<Scalar> res;
2049 diag_expr(_expr<E>
const& u) : _u(u)
2051 GISMO_ASSERT(0== _u.cols()%_u.rows(),
"Expecting square-block expression, got "
2052 << _u.rows() <<
" x "<< _u.cols() );
2055 const gsMatrix<Scalar> & eval(
const index_t k)
const
2060 const index_t r = _u.cols() / cb;
2062 for (
index_t i = 0; i!=r; ++i)
2063 res.row(i) = tmp.middleCols(i*cb,cb).diagonal();
2067 index_t rows()
const {
return _u.cols() / _u.rows(); }
2068 index_t cols()
const {
return _u.rows(); }
2070 void parse(gsExprHelper<Scalar> & evList)
const
2071 { _u.parse(evList); }
2073 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2074 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
2077 void print(std::ostream &os)
const { os <<
"diag("; _u.print(os); os<<
")"; }
2081 template <
typename E> EIGEN_STRONG_INLINE
2083 {
return diag_expr<E>(u); }
2085 #define GISMO_EXPR_VECTOR_EXPRESSION(name, mname, isSv) \
2086 template<class E> class name##_##expr : public _expr<name##_##expr<E> > { \
2087 typename E::Nested_t _u; \
2089 typedef typename E::Scalar Scalar; \
2090 enum {Space= E::Space, ScalarValued= isSv, ColBlocks= E::ColBlocks}; \
2091 name##_##expr(_expr<E> const& u) : _u(u) { } \
2092 mutable Temporary_t tmp; \
2093 const Temporary_t & eval(const index_t k) const { \
2094 tmp = _u.eval(k).mname(); return tmp; } \
2095 index_t rows() const { return isSv ? 0 : _u.rows(); } \
2096 index_t cols() const { return isSv ? 0 : _u.cols(); } \
2097 void parse(gsExprHelper<Scalar> & evList) const { _u.parse(evList); } \
2098 const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();} \
2099 const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();} \
2100 void print(std::ostream &os) const \
2101 { os << #name <<"("; _u.print(os); os <<")"; } \
2122 #undef GISMO_EXPR_VECTOR_EXPRESSION
2128 class asdiag_expr :
public _expr<asdiag_expr<E> >
2131 typedef typename E::Scalar Scalar;
2133 typename E::Nested_t _u;
2137 enum{Space = E::Space, ScalarValued= 0, ColBlocks= E::ColBlocks};
2139 asdiag_expr(_expr<E>
const& u) : _u(u) { }
2143 const gsMatrix<Scalar> & eval(
const index_t k)
const
2145 auto m = _u.eval(k);
2149 for (
index_t i = 0; i!=c; ++i)
2150 res.middleCols(i*r,r) = m.col(i).asDiagonal();
2154 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2155 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
2157 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
2159 index_t rows()
const {
return _u.rows(); }
2160 index_t cols()
const {
return _u.rows() * _u.cols(); }
2161 void parse(gsExprHelper<Scalar> & evList)
const
2162 { _u.parse(evList); }
2164 void print(std::ostream &os)
const { os <<
"diag("; _u.print(os); os <<
")";}
2169 class max_expr :
public _expr<max_expr<E> >
2172 typedef typename E::Scalar Scalar;
2173 enum {ScalarValued = 0, Space = E::Space, ColBlocks = 1};
2175 typename E::Nested_t _u;
2176 mutable gsMatrix<Scalar> tmp;
2177 mutable gsMatrix<Scalar> res;
2181 max_expr(_expr<E>
const& u) : _u(u)
2186 const gsMatrix<Scalar> & eval(
const index_t k)
const {
return eval_impl(_u,k); }
2188 index_t rows()
const {
return 1; }
2189 index_t cols()
const {
return 1; }
2190 void setFlag()
const { _u.setFlag(); }
2192 void parse(gsExprHelper<Scalar> & evList)
const
2193 { _u.parse(evList); }
2195 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2196 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
2197 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
2199 void print(std::ostream &os)
const { os <<
"max("; _u.print(os); os<<
")"; }
2201 template<
class U>
inline
2202 typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value,
const gsMatrix<Scalar> & >::type
2203 eval_impl(
const U & u,
const index_t k)
const
2207 res.resize(1,u.cardinality());
2209 for (
index_t c=0; c!=_u.cardinality(); c++)
2210 res(0,c) = tmp.block(0,c*u.cols(),u.rows(),u.cols()).maxCoeff();
2212 for (
index_t c=0; c!=_u.rows(); c++)
2213 res(0,c) = tmp.block(c*u.rows(),0,u.rows(),u.cols()).maxCoeff();
2218 template<
class U>
inline
2219 typename util::enable_if< !util::is_same<U,gsFeSpace<Scalar> >::value,
const gsMatrix<Scalar> & >::type
2220 eval_impl(
const U & u,
const index_t k)
const
2222 res = u.eval(k).colwise().maxCoeff();
2228 class rowsum_expr :
public _expr<rowsum_expr<E> >
2231 typedef typename E::Scalar Scalar;
2232 enum {ScalarValued = 0, Space = E::Space, ColBlocks = E::ColBlocks};
2234 typename E::Nested_t _u;
2235 mutable gsMatrix<Scalar> tmp;
2239 rowsum_expr(_expr<E>
const& u) : _u(u)
2244 const gsMatrix<Scalar> & eval(
const index_t k)
const
2246 tmp = _u.eval(k).rowwise().sum();
2250 index_t rows()
const {
return _u.rows(); }
2251 index_t cols()
const {
return 1; }
2252 void setFlag()
const { _u.setFlag(); }
2254 void parse(gsExprHelper<Scalar> & evList)
const
2255 { _u.parse(evList); }
2257 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2258 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
2259 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
2261 void print(std::ostream &os)
const { os <<
"rowsum("; _u.print(os); os<<
")"; }
2265 class colsum_expr :
public _expr<colsum_expr<E> >
2268 typedef typename E::Scalar Scalar;
2269 enum {ScalarValued = 0, Space = E::Space, ColBlocks = E::ColBlocks};
2271 typename E::Nested_t _u;
2272 mutable gsMatrix<Scalar> tmp;
2276 colsum_expr(_expr<E>
const& u) : _u(u)
2281 const gsMatrix<Scalar> & eval(
const index_t k)
const
2283 tmp = _u.eval(k).colwise().sum();
2287 index_t rows()
const {
return _u.rows(); }
2288 index_t cols()
const {
return 1; }
2289 void setFlag()
const { _u.setFlag(); }
2291 void parse(gsExprHelper<Scalar> & evList)
const
2292 { _u.parse(evList); }
2294 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2295 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
2296 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
2298 void print(std::ostream &os)
const { os <<
"colsum("; _u.print(os); os<<
")"; }
2307 typedef real_t Scalar;
2308 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
2322 index_t rows()
const {
return _dim; }
2323 index_t cols()
const {
return _dim; }
2329 void print(std::ostream &os)
const { os <<
"id("<<_dim <<
")";}
2338 typedef real_t Scalar;
2339 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
2353 index_t rows()
const {
return _mat.rows(); }
2354 index_t cols()
const {
return _mat.cols(); }
2360 void print(std::ostream &os)
const { os <<
"constMat";}
2367 class sign_expr :
public _expr<sign_expr<E> >
2369 typename E::Nested_t _u;
2370 typename E::Scalar _tol;
2372 typedef typename E::Scalar Scalar;
2373 enum {ScalarValued = 1, Space = E::Space, ColBlocks= 0};
2375 sign_expr(_expr<E>
const& u, Scalar tolerance = 0.0) : _u(u),_tol(tolerance){
2376 GISMO_ASSERT( _tol >= 0,
"Tolerance for sign_expr should be a positive number.");
2379 Scalar eval(
const index_t k)
const
2381 const Scalar v = _u.val().eval(k);
2382 return ( v>_tol ? 1 : ( v<-_tol ? -1 : 0 ) );
2385 static index_t rows() {
return 0; }
2386 static index_t cols() {
return 0; }
2388 void parse(gsExprHelper<Scalar> & el)
const
2391 static bool isScalar() {
return true; }
2393 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
2394 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
2396 void print(std::ostream &os)
const { os<<
"sgn("; _u.print(os); os <<
")"; }
2403 class exp_expr :
public _expr<exp_expr<E> >
2405 typename E::Nested_t _u;
2407 typedef typename E::Scalar Scalar;
2408 enum {ScalarValued = 1, Space = E::Space, ColBlocks= 0};
2410 exp_expr(_expr<E>
const& u) : _u(u) { }
2412 Scalar eval(
const index_t k)
const
2414 const Scalar v = _u.val().eval(k);
2415 return math::exp(v);
2418 static index_t rows() {
return 0; }
2419 static index_t cols() {
return 0; }
2421 void parse(gsExprHelper<Scalar> & el)
const
2424 static bool isScalar() {
return true; }
2426 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
2427 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
2429 void print(std::ostream &os)
const { os<<
"exp("; _u.print(os); os <<
")"; }
2436 class ppart_expr :
public _expr<ppart_expr<E> >
2439 typedef typename E::Scalar Scalar;
2440 enum {ScalarValued = E::ScalarValued, Space = E::Space, ColBlocks= E::ColBlocks};
2442 typename E::Nested_t _u;
2443 mutable gsMatrix<Scalar> res;
2446 ppart_expr(_expr<E>
const& u) : _u(u) { }
2448 const gsMatrix<Scalar> & eval(
index_t k)
const
2450 res = _u.eval(k).cwiseMax(0.0);
2455 index_t rows()
const {
return _u.rows(); }
2456 index_t cols()
const {
return _u.cols(); }
2458 void parse(gsExprHelper<Scalar> & el)
const
2461 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar();}
2462 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar();}
2464 void print(std::ostream &os)
const { os<<
"posPart("; _u.print(os); os <<
")"; }
2471 class ppartval_expr :
public _expr<ppartval_expr<E> >
2473 typename E::Nested_t _u;
2475 typedef typename E::Scalar Scalar;
2476 enum {ScalarValued = 1, Space = 0, ColBlocks= 0};
2480 ppartval_expr(_expr<E>
const& u) : _u(u) { }
2482 Scalar & eval(
index_t k)
const
2484 res = std::max(0.0,_u.eval(k));
2488 index_t rows()
const {
return 0; }
2489 index_t cols()
const {
return 0; }
2491 void parse(gsExprHelper<Scalar> & evList)
const
2492 { _u.parse(evList); }
2494 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
2495 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
2497 void print(std::ostream &os)
const { os<<
"posPart("; _u.print(os); os <<
")"; }
2504 class pow_expr :
public _expr<pow_expr<E> >
2506 typename E::Nested_t _u;
2509 typedef typename E::Scalar Scalar;
2510 enum {ScalarValued = 1, Space = E::Space, ColBlocks= E::ColBlocks};
2514 pow_expr(_expr<E>
const& u, Scalar q) : _u(u), _q(q) { }
2516 Scalar eval(
const index_t k)
const
2518 const Scalar v = _u.val().eval(k);
2519 return math::pow(v,_q);
2522 static index_t rows() {
return 0; }
2523 static index_t cols() {
return 0; }
2525 void parse(gsExprHelper<Scalar> & el)
const
2528 static bool isScalar() {
return true; }
2530 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
2531 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
2533 void print(std::ostream &os)
const { os<<
"pow("; _u.print(os); os <<
")"; }
2541 template <
typename E1,
typename E2>
2545 typedef typename E1::Scalar Scalar;
2546 enum {ScalarValued = 0, ColBlocks = 1};
2547 enum {Space = E2::Space};
2550 typename E1::Nested_t _u;
2551 typename E2::Nested_t _v;
2562 const index_t N = _v.cols() / (r*r);
2564 const auto uEv = _u.eval(k);
2565 const auto vEv = _v.eval(k);
2567 res.resize(r, N*r*r);
2569 for (
index_t s = 0; s!=r; ++s)
2570 for (
index_t i = 0; i!=N; ++i)
2572 res.middleCols((s*N + i)*r,r).noalias() =
2573 uEv.col(s) * vEv.middleCols((s*N + i)*r,r).row(s);
2580 index_t rows()
const {
return _u.cols(); }
2581 index_t cols()
const {
return _v.cols(); }
2584 { _u.parse(evList); _v.parse(evList); }
2589 void print(std::ostream &os)
const { os <<
"matrix_by_space("; _u.print(os); os<<
")"; }
2597 template <
typename E1,
typename E2>
2601 typedef typename E1::Scalar Scalar;
2602 enum {ScalarValued = 0, ColBlocks = 1};
2603 enum {Space = E2::Space};
2606 typename E1::Nested_t _u;
2607 typename E2::Nested_t _v;
2618 const index_t N = _v.cols() / (r*r);
2620 const auto uEv = _u.eval(k);
2621 const auto vEv = _v.eval(k);
2623 res.resize(r, N*r*r);
2624 for (
index_t s = 0; s!=r; ++s)
2625 for (
index_t i = 0; i!=N; ++i)
2627 res.middleCols((s*N + i)*r,r).noalias() =
2628 uEv.transpose()*vEv.middleCols((s*N + i)*r,r).transpose();
2634 index_t rows()
const {
return _u.cols(); }
2635 index_t cols()
const {
return _v.cols(); }
2638 { _u.parse(evList); _v.parse(evList); }
2643 void print(std::ostream &os)
const { os <<
"matrix_by_space_tr("; _u.print(os); os<<
")"; }
2650 class value_expr :
public _expr<value_expr<E> >
2652 typename E::Nested_t _u;
2655 typedef typename E::Scalar Scalar;
2656 value_expr(_expr<E>
const& u) : _u(u)
2663 enum {Space= 0, ScalarValued= 1, ColBlocks= 0};
2665 Scalar eval(
const index_t k)
const {
return eval_impl(_u,k); }
2668 inline value_expr<E> val()
const {
return *
this; }
2669 index_t rows()
const {
return 0; }
2670 index_t cols()
const {
return 0; }
2671 void parse(gsExprHelper<Scalar> & evList)
const
2672 { _u.parse(evList); }
2674 static bool isScalar() {
return true; }
2676 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
2677 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
2679 void print(std::ostream &os)
const { _u.print(os); }
2684 template<
class U>
static inline
2685 typename util::enable_if<U::ScalarValued,Scalar>::type
2686 eval_impl(
const U & u,
const index_t k) {
return u.eval(k); }
2688 template<
class U>
static inline
2689 typename util::enable_if<!U::ScalarValued,Scalar>::type
2690 eval_impl(
const U & u,
const index_t k) {
return u.eval(k).value(); }
2694 class abs_expr :
public _expr<abs_expr<E> >
2696 typename E::Nested_t _u;
2699 typedef typename E::Scalar Scalar;
2700 explicit abs_expr(_expr<E>
const& u) : _u(u) { }
2703 enum {Space= 0, ScalarValued= 1, ColBlocks= 0};
2705 Scalar eval(
const index_t k)
const {
return abs_expr::eval_impl(_u,k); }
2707 index_t rows()
const {
return _u.rows(); }
2708 index_t cols()
const {
return _u.cols(); }
2709 void parse(gsExprHelper<Scalar> & evList)
const
2710 { _u.parse(evList); }
2712 static bool isScalar() {
return true; }
2714 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get();}
2715 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get();}
2717 void print(std::ostream &os)
const { _u.print(os); }
2722 template<
class U>
static inline
2723 typename util::enable_if<U::ScalarValued,Scalar>::type
2725 template<
class U>
static inline
2726 typename util::enable_if<!U::ScalarValued,gsMatrix<Scalar> >::type
2727 eval_impl(
const U & u,
const index_t k) {
return u.eval(k).cwiseAbs(); }
2736 class grad_expr :
public _expr<grad_expr<E> >
2738 typename E::Nested_t _u;
2740 enum {Space = E::Space, ScalarValued= 0, ColBlocks= 0};
2742 typedef typename E::Scalar Scalar;
2743 mutable gsMatrix<Scalar> tmp;
2745 grad_expr(
const E & u) : _u(u)
2746 {
GISMO_ASSERT(1==u.dim(),
"grad(.) requires 1D variable, use jac(.) instead.");}
2748 const gsMatrix<Scalar> & eval(
const index_t k)
const
2753 tmp = _u.data().values[1].reshapeCol(k, cols(), cardinality_impl()).transpose();
2757 index_t rows()
const {
return 1 ; }
2759 index_t cols()
const {
return _u.source().domainDim(); }
2761 index_t cardinality_impl()
const
2762 {
return _u.data().values[1].rows() / cols(); }
2764 void parse(gsExprHelper<Scalar> & evList)
const
2767 _u.data().flags |= NEED_GRAD;
2770 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2771 const gsFeSpace<Scalar> & colVar()
const
2772 {
return gsNullExpr<Scalar>::get();}
2774 void print(std::ostream &os)
const { os <<
"\u2207("; _u.print(os); os <<
")"; }
2777 template<
class U>
static inline
2778 typename util::enable_if<util::is_same<U,gsComposition<Scalar> >::value,
2779 const gsMatrix<Scalar> &>::type
2780 eval_impl(
const U & u,
const index_t k)
2794 class grad_expr<gsFeSolution<T> > :
public _expr<grad_expr<gsFeSolution<T> > >
2797 const gsFeSolution<T> _u;
2801 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
2803 explicit grad_expr(
const gsFeSolution<T> & u) : _u(u) { }
2805 mutable gsMatrix<T> res;
2806 const gsMatrix<T> & eval(
index_t k)
const
2808 GISMO_ASSERT(1==_u.data().actives.cols(),
"Single actives expected");
2810 res.setZero(_u.dim(), _u.parDim());
2811 const gsDofMapper & map = _u.mapper();
2812 for (
index_t c = 0; c!= _u.dim(); c++)
2814 for (
index_t i = 0; i!=_u.data().actives.size(); ++i)
2816 const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId,c);
2817 if ( map.is_free_index(ii) )
2819 res.row(c) += _u.coefs().at(ii) *
2820 _u.data().values[1].col(k).segment(i*_u.parDim(), _u.parDim()).transpose();
2825 _u.fixedPart().at( map.global_to_bindex(ii) ) *
2826 _u.data().values[1].col(k).segment(i*_u.parDim(), _u.parDim()).transpose();
2833 index_t rows()
const {
return _u.dim();}
2834 index_t cols()
const {
return _u.parDim(); }
2836 const gsFeSpace<Scalar> & rowVar()
const
2837 {
return gsNullExpr<Scalar>::get();}
2838 const gsFeSpace<Scalar> & colVar()
const
2839 {
return gsNullExpr<Scalar>::get();}
2841 void parse(gsExprHelper<Scalar> & evList)
const
2844 evList.add(_u.space());
2848 void print(std::ostream &os)
const { os <<
"\u2207(s)"; }
2858 class dJacdc_expr :
public _expr<dJacdc_expr<E> >
2860 typename E::Nested_t _u;
2862 enum{ Space = E::Space, ScalarValued = 0, ColBlocks = (1==E::Space?1:0)};
2864 typedef typename E::Scalar Scalar;
2866 mutable gsMatrix<Scalar> res;
2869 dJacdc_expr(
const E & u,
index_t c) : _u(u), _c(c)
2870 {
GISMO_ASSERT(1==u.dim(),
"grad(.) requires 1D variable, use jac(.) instead.");}
2872 const gsMatrix<Scalar> & eval(
const index_t k)
const
2874 index_t dd = _u.source().domainDim();
2876 res.setZero(dd, dd*n);
2878 gsMatrix<Scalar>
grad = _u.data().values[1].reshapeCol(k, dd, n);
2879 for(
index_t i = 0; i < n; i++){
2880 res.row(_c).segment(i*dd,dd) = grad.col(i);
2885 index_t rows()
const {
return _u.source().domainDim(); }
2887 index_t cols()
const {
return _u.source().domainDim()*_u.rows(); }
2889 void parse(gsExprHelper<Scalar> & evList)
const
2892 _u.data().flags |= NEED_GRAD;
2895 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
2896 const gsFeSpace<Scalar> & colVar()
const
2897 {
return gsNullExpr<Scalar>::get();}
2899 void print(std::ostream &os)
const { os <<
"dJacdc("; _u.print(os); os <<
")"; }
2907 class nabla_expr :
public _expr<nabla_expr<T> >
2909 typename gsFeVariable<T>::Nested_t u;
2920 nabla_expr(
const gsFeVariable<T> & _u) : u(_u)
2926 mutable gsMatrix<Scalar> res;
2928 const gsMatrix<Scalar> & eval(
const index_t k)
const
2932 res.setZero(rows(), d);
2934 for (
index_t i = 0; i!=d; ++i)
2935 res.col(i).segment(i*n,n) = u.data().values[1].reshapeCol(k, d, n).row(i);
2939 index_t rows()
const {
return u.rows(); }
2940 index_t cols()
const {
return u.cols(); }
2942 void parse(gsExprHelper<Scalar> & evList)
const
2945 u.data().flags |= NEED_GRAD;
2948 const gsFeSpace<T> & rowVar()
const {
return u.rowVar(); }
2949 const gsFeSpace<T> & colVar()
const
2950 {
return gsNullExpr<Scalar>::get();}
2952 void print(std::ostream &os)
const { os <<
"nabla("; u.print(os); os <<
")"; }
2962 class nabla2_expr :
public _expr<nabla2_expr<T> >
2964 typename gsFeVariable<T>::Nested_t u;
2975 nabla2_expr(
const gsFeVariable<T> & _u)
2982 return u.data().values[2]
2983 .reShapeCol(k, u.data().values[2].rows()/u.cSize(), u.cSize() )
2984 .topRows(u.parDim()).transpose();
2987 index_t rows()
const {
return u.rows(); }
2988 index_t cols()
const {
return u.parDim(); }
2990 void parse(gsExprHelper<Scalar> & evList)
const
2996 const gsFeSpace<T> & rowVar()
const {
return u.rowVar(); }
2997 const gsFeSpace<T> & colVar()
const
2998 {
return gsNullExpr<T>::get();}
3013 typename gsGeometryMap<T>::Nested_t _G;
3017 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3019 explicit onormal_expr(
const gsGeometryMap<T> & G) : _G(G) { }
3021 auto eval(
const index_t k)
const -> decltype(_G.data().outNormals.col(k))
3022 {
return _G.data().outNormals.col(k); }
3024 index_t rows()
const {
return _G.source().targetDim(); }
3025 index_t cols()
const {
return 1; }
3027 const gsFeSpace<T> & rowVar()
const {
return gsNullExpr<T>::get();}
3028 const gsFeSpace<T> & colVar()
const {
return gsNullExpr<T>::get();}
3036 void print(std::ostream &os)
const { os <<
"nv("; _G.print(os); os <<
")"; }
3046 typename gsGeometryMap<T>::Nested_t _G;
3050 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3054 GISMO_ENSURE( _G.source().domainDim()+1 == _G.source().targetDim(),
"Surface normal requires codimension 1");
3057 auto eval(
const index_t k)
const -> decltype(_G.data().normals.col(k))
3058 {
return _G.data().normals.col(k); }
3060 index_t rows()
const {
return _G.source().targetDim(); }
3061 index_t cols()
const {
return 1; }
3063 const gsFeSpace<T> & rowVar()
const {
return gsNullExpr<T>::get();}
3064 const gsFeSpace<T> & colVar()
const {
return gsNullExpr<T>::get();}
3072 void print(std::ostream &os)
const { os <<
"sn("; _G.print(os); os <<
")"; }
3082 typename gsGeometryMap<T>::Nested_t _G;
3086 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3093 if (_G.targetDim()==2)
3095 res = _G.data().outNormals.col(k);
3096 std::swap( res(0,0), res(1,0) );
3100 else if (_G.targetDim()==3)
3103 res.col3d(0) = _G.data().normals.col3d(k)
3104 .cross( _G.data().outNormals.col3d(k) );
3108 GISMO_ERROR(
"Function not implemented for dimension"<<_G.targetDim());
3112 index_t rows()
const {
return _G.source().targetDim(); }
3113 index_t cols()
const {
return 1; }
3125 void print(std::ostream &os)
const { os <<
"tv("; _G.print(os); os <<
")"; }
3134 typename E::Nested_t _u;
3137 typedef typename E::Scalar Scalar;
3138 enum {Space = E::Space, ScalarValued= 0, ColBlocks= 0};
3142 auto eval(
const index_t k)
const -> decltype(_u.data().laplacians.col(k))
3145 return _u.data().laplacians.col(k);
3151 index_t rows()
const {
return _u.data().laplacians.rows(); }
3152 index_t cols()
const {
return 1; }
3154 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
3165 void print(std::ostream &os)
const { os <<
"\u2206("; _u.print(os); os <<
")"; }
3172 class lapl_expr<gsFeSolution<T> > :
public _expr<lapl_expr<gsFeSolution<T> > >
3175 const gsFeSolution<T> _u;
3179 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3181 lapl_expr(
const gsFeSolution<T> & u) : _u(u) { }
3183 mutable gsMatrix<T> res;
3184 const gsMatrix<T> & eval(
const index_t k)
const
3186 GISMO_ASSERT(1==_u.data().actives.cols(),
"Single actives expected");
3188 res.setZero(_u.dim(), 1);
3189 const gsDofMapper & map = _u.mapper();
3191 index_t numActs = _u.data().values[0].rows();
3192 index_t numDers = _u.parDim() * (_u.parDim() + 1) / 2;
3195 for (
index_t c = 0; c!= _u.dim(); c++)
3196 for (
index_t i = 0; i!=numActs; ++i)
3198 const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId,c);
3199 deriv2 = _u.data().values[2].block(i*numDers,k,_u.parDim(),1);
3200 if ( map.is_free_index(ii) )
3201 res.at(c) += _u.coefs().at(ii) * deriv2.sum();
3203 res.at(c) +=_u.fixedPart().at( map.global_to_bindex(ii) ) * deriv2.sum();
3208 index_t rows()
const {
return _u.dim(); }
3209 index_t cols()
const {
return 1; }
3211 void parse(gsExprHelper<Scalar> & evList)
const
3213 evList.add(_u.space());
3217 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<T>::get();}
3218 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<T>::get();}
3220 void print(std::ostream &os)
const { os <<
"\u2206(s)"; }
3227 class fform2nd_expr :
public _expr<fform2nd_expr<T> >
3229 typename gsGeometryMap<T>::Nested_t _G;
3232 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
3234 fform2nd_expr(
const gsGeometryMap<T> & G) : _G(G) { }
3236 const gsAsConstMatrix<Scalar> eval(
const index_t k)
const
3238 return gsAsConstMatrix<Scalar>(_G.data().fundForms.col(k).data(),rows(),cols());
3241 index_t rows()
const {
return _G.source().domainDim() ; }
3242 index_t cols()
const {
return _G.source().domainDim() ; }
3244 void parse(gsExprHelper<Scalar> & evList)
const
3250 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<T>::get();}
3251 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<T>::get();}
3253 void print(std::ostream &os)
const { os <<
"fform2nd("; _G.print(os); os <<
")"; }
3261 class jacInv_expr :
public _expr<jacInv_expr<T> >
3263 typename gsGeometryMap<T>::Nested_t _G;
3266 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
3268 jacInv_expr(
const gsGeometryMap<T> & G) : _G(G)
3274 MatExprType eval(
const index_t k)
const {
return _G.data().jacInvTr.reshapeCol(k,cols(),rows()).transpose(); }
3276 index_t rows()
const {
return _G.source().domainDim(); }
3277 index_t cols()
const {
return _G.source().targetDim(); }
3279 void parse(gsExprHelper<Scalar> & evList)
const
3285 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<T>::get();}
3286 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<T>::get();}
3291 void print(std::ostream &os)
const { os <<
"jacInv("; _G.print(os); os <<
")"; }
3298 class jac_expr :
public _expr<jac_expr<E> >
3300 typename E::Nested_t _u;
3302 enum {ColBlocks = (1==E::Space?1:0) };
3303 enum {Space = E::Space, ScalarValued= 0 };
3305 typedef typename E::Scalar Scalar;
3307 mutable gsMatrix<Scalar> res;
3309 jac_expr(
const E & _u_)
3317 res = _u.data().values[1].col(k).transpose().blockDiag(_u.dim());
3321 res = _u.data().values[1]
3322 .reshapeCol(k, _u.parDim(), _u.targetDim()).transpose()
3323 .blockDiag(_u.dim());
3328 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
3330 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get(); }
3332 index_t rows()
const {
return rows_impl(_u); }
3333 index_t cols()
const {
return cols_impl(_u); }
3338 index_t cardinality_impl()
const
3340 return _u.dim() * _u.data().actives.rows();
3343 void parse(gsExprHelper<Scalar> & evList)
const
3350 void print(std::ostream &os)
const { os <<
"\u2207("; _u.print(os);os <<
")"; }
3358 template<
class U>
inline
3359 typename util::enable_if<!(util::is_same<U,gsFeSpace<Scalar> >::value),
index_t >::type
3360 rows_impl(
const U & u)
const
3362 return u.source().targetDim();
3365 template<
class U>
inline
3366 typename util::enable_if< (util::is_same<U,gsFeSpace<Scalar> >::value),
index_t >::type
3367 rows_impl(
const U & u)
const
3372 template<
class U>
inline
3373 typename util::enable_if<!(util::is_same<U,gsFeSpace<Scalar> >::value),
index_t >::type
3374 cols_impl(
const U & u)
const
3376 return u.source().domainDim();
3379 template<
class U>
inline
3380 typename util::enable_if< (util::is_same<U,gsFeSpace<Scalar> >::value),
index_t >::type
3381 cols_impl(
const U & u)
const
3383 return u.source().domainDim();
3390 template<
class U>
inline
3391 typename util::enable_if<!(util::is_same<U,gsFeSpace<Scalar> >::value),
const gsFeSpace<Scalar> & >::type
3394 return gsNullExpr<Scalar>::get();
3397 template<
class U>
inline
3398 typename util::enable_if<(util::is_same<U,gsFeSpace<Scalar> >::value),
const gsFeSpace<Scalar> & >::type
3409 class jac_expr<gsGeometryMap<T> > :
public _expr<jac_expr<gsGeometryMap<T> > >
3411 typename gsGeometryMap<T>::Nested_t _G;
3415 enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3417 jac_expr(
const gsGeometryMap<T> & G) : _G(G) { }
3421 return _G.data().values[1]
3422 .reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
3425 index_t rows()
const {
return _G.source().targetDim(); }
3427 index_t cols()
const {
return _G.source().domainDim(); }
3429 static const gsFeSpace<Scalar> & rowVar() {
return gsNullExpr<Scalar>::get(); }
3430 static const gsFeSpace<Scalar> & colVar() {
return gsNullExpr<Scalar>::get(); }
3432 void parse(gsExprHelper<Scalar> & evList)
const
3438 meas_expr<T> absDet()
const
3440 GISMO_ASSERT(rows() == cols(),
"The Jacobian matrix is not square");
3441 return meas_expr<T>(_G);
3444 jacInv_expr<T> inv()
const
3446 GISMO_ASSERT(rows() == cols(),
"The Jacobian matrix is not square");
3447 return jacInv_expr<T>(_G);
3451 jacInv_expr<T> ginv()
const {
return jacInv_expr<T>(_G); }
3453 void print(std::ostream &os)
const { os <<
"\u2207("; _G.print(os); os <<
")"; }
3457 class hess_expr :
public _expr<hess_expr<E> >
3460 typedef typename E::Scalar Scalar;
3462 typename E::Nested_t _u;
3463 mutable gsMatrix<Scalar> res;
3465 enum {ScalarValued = 0, ColBlocks = (1==E::Space?1:0) };
3466 enum {Space = E::Space };
3469 hess_expr(
const E & u) : _u(u)
3475 const gsMatrix<Scalar> & eval(
const index_t k)
const
3477 const gsFuncData<Scalar> & dd = _u.data();
3478 const index_t sz = cardinality_impl();
3479 res.resize(dd.dim.first, sz*dd.dim.first);
3480 secDerToHessian(dd.values[2].col(k), dd.dim.first, res);
3481 res.resize(dd.dim.first, res.cols()*dd.dim.first);
3487 index_t rows()
const {
return _u.data().dim.first; }
3493 index_t cardinality_impl()
const
3495 return 2*_u.data().values[2].rows()/
3496 (_u.data().dim.first*(1+_u.data().dim.first));
3500 void parse(gsExprHelper<Scalar> & evList)
const
3503 _u.data().flags |= NEED_2ND_DER;
3506 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
3507 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get(); }
3509 void print(std::ostream &os)
const
3511 { os <<
"\u210D(U)"; }
3515 class hess_expr<gsFeSolution<T> > :
public _expr<hess_expr<gsFeSolution<T> > >
3518 const gsFeSolution<T> _u;
3522 enum{Space = 0, ScalarValued = 0, ColBlocks = 0 };
3524 hess_expr(
const gsFeSolution<T> & u) : _u(u) { }
3526 mutable gsMatrix<T> res;
3527 const gsMatrix<T> & eval(
const index_t k)
const
3529 GISMO_ASSERT(1==_u.data().actives.cols(),
"Single actives expected. Actives: \n"<<_u.data().actives);
3531 const gsDofMapper & map = _u.mapper();
3533 const index_t numActs = _u.data().values[0].rows();
3534 const index_t pdim = _u.parDim();
3535 index_t numDers = pdim*(pdim+1)/2;
3541 res.setZero(numDers,1);
3542 for (
index_t i = 0; i!=numActs; ++i)
3544 const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId,0);
3545 deriv2 = _u.data().values[2].block(i*numDers,k,numDers,1);
3546 if ( map.is_free_index(ii) )
3547 res += _u.coefs().at(ii) * deriv2;
3549 res +=_u.fixedPart().at( map.global_to_bindex(ii) ) * deriv2;
3551 secDerToHessian(res, pdim, deriv2);
3553 res.resize(pdim,pdim);
3558 res.setZero(rows(), numDers);
3559 for (
index_t c = 0; c != _u.dim(); c++)
3560 for (
index_t i = 0; i != numActs; ++i) {
3561 const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId, c);
3562 deriv2 = _u.space().data().values[2].block(i * numDers, k, numDers,
3564 if (map.is_free_index(ii))
3565 res.row(c) += _u.coefs().at(ii) * deriv2;
3567 res.row(c) += _u.fixedPart().at(map.global_to_bindex(ii)) * deriv2;
3586 return _u.parDim() * (_u.parDim() + 1) / 2;
3589 const gsFeSpace<Scalar> & rowVar()
const {
return gsNullExpr<Scalar>::get(); }
3590 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get(); }
3592 void parse(gsExprHelper<Scalar> & evList)
const
3595 evList.add(_u.space());
3599 void print(std::ostream &os)
const { os <<
"\u210D(s)"; }
3608 class dJacG_expr :
public _expr<dJacG_expr<T> >
3610 typename gsGeometryMap<T>::Nested_t _G;
3612 mutable gsMatrix<T> res;
3616 dJacG_expr(
const gsGeometryMap<T> & G) : _G(G) { }
3620 const index_t sz = _G.data().values[0].rows();
3621 const index_t s = _G.data().derivSize();
3623 res.resize(_G.data().dim.second, sz*_G.data().dim.first);
3628 index_t rows()
const {
return _G.source().targetDim(); }
3629 index_t cols()
const {
return _G.source().domainDim(); }
3631 void parse(gsExprHelper<Scalar> & evList)
const
3634 _G.data().flags |= NEED_2ND_DER;
3643 class curl_expr :
public _expr<curl_expr<T> >
3648 typename gsFeVariable<T>::Nested_t _u;
3649 mutable gsMatrix<Scalar> res;
3651 enum{ Space = 1, ScalarValued= 0, ColBlocks= 0};
3653 curl_expr(
const gsFeVariable<T> & u) : _u(u)
3654 {
GISMO_ASSERT(3==u.dim(),
"curl(.) requires 3D variable."); }
3656 const gsMatrix<T> & eval(
const index_t k)
const
3658 res.setZero( rows(), _u.dim());
3659 const index_t na = _u.data().values[0].rows();
3660 gsAsConstMatrix<T, Dynamic, Dynamic> pd =
3661 _u.data().values[1].reshapeCol(k, cols(), na);
3663 res.col(0).segment(na ,na) = -pd.row(2);
3664 res.col(0).segment(2*na,na) = pd.row(1);
3665 res.col(1).segment(0 ,na) = pd.row(2);
3666 res.col(1).segment(2*na,na) = -pd.row(0);
3667 res.col(2).segment(0 ,na) = -pd.row(1);
3668 res.col(2).segment(na ,na) = pd.row(0);
3672 index_t rows()
const {
return _u.dim() * _u.data().values[0].rows(); }
3673 index_t cols()
const {
return _u.data().dim.first; }
3675 void parse(gsExprHelper<Scalar> & evList)
const
3678 _u.data().flags |= NEED_GRAD;
3681 const gsFeSpace<T> & rowVar()
const {
return _u.rowVar(); }
3682 const gsFeSpace<T> & colVar()
const {
return gsNullExpr<T>::get();}
3684 void print(std::ostream &os)
const { os <<
"curl("; _u.print(os); os <<
")"; }
3697 template <
typename E1,
typename E2>
3698 class mult_expr<E1,E2,false> :
public _expr<mult_expr<E1, E2, false> >
3700 typename E1::Nested_t _u;
3701 typename E2::Nested_t _v;
3704 enum {ScalarValued = E1::ScalarValued && E2::ScalarValued,
3705 ColBlocks = E2::ColBlocks};
3706 enum {Space = (int)E1::Space + (
int)E2::Space };
3708 typedef typename E1::Scalar Scalar;
3710 mult_expr(_expr<E1>
const& u,
3714 mutable Temporary_t tmp;
3715 const Temporary_t & eval(
const index_t k)
const
3717 GISMO_ASSERT(0==_u.cols()*_v.rows() || _u.cols() == _v.rows(),
3718 "Wrong dimensions "<<_u.cols()<<
"!="<<_v.rows()<<
" in * operation:\n"
3719 << _u <<
" times \n" << _v );
3722 tmp = _u.eval(k) * _v.eval(k);
3726 index_t rows()
const {
return E1::ScalarValued ? _v.rows() : _u.rows(); }
3727 index_t cols()
const {
return E2::ScalarValued ? _u.cols() : _v.cols(); }
3728 void parse(gsExprHelper<Scalar> & evList)
const
3729 { _u.parse(evList); _v.parse(evList); }
3732 index_t cardinality_impl()
const
3733 {
return 0==E1::Space ? _v.cardinality(): _u.cardinality(); }
3735 const gsFeSpace<Scalar> & rowVar()
const
3736 {
return 0==E1::Space ? _v.rowVar() : _u.rowVar(); }
3737 const gsFeSpace<Scalar> & colVar()
const
3738 {
return 0==E2::Space ? _u.colVar() : _v.colVar(); }
3740 void print(std::ostream &os)
const { _u.print(os); os<<
"*"; _v.print(os); }
3757 template <
typename E1,
typename E2>
3758 class mult_expr<E1, E2, true> :
public _expr<mult_expr<E1, E2, true> >
3761 typedef typename E2::Scalar Scalar;
3763 typename E1::Nested_t _u;
3764 typename E2::Nested_t _v;
3766 mutable gsMatrix<Scalar> res;
3768 enum {ScalarValued = 0, ColBlocks = E1::ColBlocks};
3769 enum {Space = (int)E1::Space + (
int)E2::Space };
3771 mult_expr(_expr<E1>
const& u,
3778 const gsMatrix<Scalar> & eval(
const index_t k)
const
3782 const index_t nb = _u.cardinality();
3783 const auto tmpA = _u.eval(k);
3784 const auto tmpB = _v.eval(k);
3789 if ( 1 == _v.cardinality() )
3791 res.resize(ur, vc*nb);
3792 GISMO_ASSERT(tmpA.cols()==uc*nb,
"Dimension error.. "<< tmpA.cols()<<
"!="<<uc*nb );
3795 for (
index_t i = 0; i!=nb; ++i)
3796 res.middleCols(i*vc,vc).noalias()
3797 = tmpA.middleCols(i*uc,uc) * tmpB;
3804 const index_t nbv = _v.cardinality();
3805 res.resize(ur*nb, vc*nbv);
3806 for (
index_t i = 0; i!=nb; ++i)
3807 for (
index_t j = 0; j!=nbv; ++j)
3809 res.block(i*ur,j*vc,ur,vc).noalias() =
3810 tmpA.middleCols(i*uc,uc) * tmpB.middleCols(j*vc,vc);
3825 void parse(gsExprHelper<Scalar> & evList)
const
3826 { _u.parse(evList); _v.parse(evList); }
3829 index_t cardinality_impl()
const {
return _u.cardinality(); }
3831 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
3832 const gsFeSpace<Scalar> & colVar()
const
3834 if ( 1 == _v.cardinality() )
3840 void print(std::ostream &os)
const
3841 { os <<
"("; _u.print(os);os <<
"*";_v.print(os);os <<
")"; }
3849 template <
typename E2>
3850 class mult_expr<typename E2::Scalar, E2, false>
3851 :
public _expr<mult_expr<typename E2::Scalar, E2, false> >
3855 typedef typename E2::Scalar Scalar;
3858 typename E2::Nested_t _v;
3862 enum {ScalarValued = E2::ScalarValued, ColBlocks = E2::ColBlocks};
3863 enum {Space = E2::Space};
3865 mult_expr(Scalar
const & c, _expr<E2>
const& v)
3868 EIGEN_STRONG_INLINE AutoReturn_t eval(
const index_t k)
const
3870 return ( _c * _v.eval(k) );
3874 index_t rows()
const {
return _v.rows(); }
3875 index_t cols()
const {
return _v.cols(); }
3877 void parse(gsExprHelper<Scalar> & evList)
const
3878 { _v.parse(evList); }
3880 index_t cardinality_impl()
const
3881 {
return _v.cardinality(); }
3883 const gsFeSpace<Scalar> & rowVar()
const {
return _v.rowVar(); }
3884 const gsFeSpace<Scalar> & colVar()
const {
return _v.colVar(); }
3886 void print(std::ostream &os)
const { os << _c <<
"*";_v.print(os); }
3890 template <
typename E1,
typename E2>
3891 class collapse_expr :
public _expr<collapse_expr<E1, E2> >
3893 typename E1::Nested_t _u;
3894 typename E2::Nested_t _v;
3897 enum {ScalarValued = 0, ColBlocks = 0};
3898 enum { Space = (int)E1::Space + (
int)E2::Space };
3900 typedef typename E1::Scalar Scalar;
3902 mutable gsMatrix<Scalar> res;
3904 collapse_expr(_expr<E1>
const& u,
3909 const gsMatrix<Scalar> &
3913 const auto tmpA = _u.eval(k);
3914 const auto tmpB = _v.eval(k);
3920 for (
index_t i = 0; i!=nb; ++i)
3922 res.row(i).transpose().noalias() = tmpA.middleCols(i*ur,ur) * tmpB;
3925 else if (E2::ColBlocks)
3929 for (
index_t i = 0; i!=nb; ++i)
3931 res.row(i).noalias() = tmpA * tmpB.middleCols(i*ur,ur);
3938 index_t rows()
const {
return E1::ColBlocks ? _u.cols() / _v.rows() : _v.cols() / _u.cols() ; }
3939 index_t cols()
const {
return E1::ColBlocks ? _v.rows() : _u.cols(); }
3941 void parse(gsExprHelper<Scalar> & evList)
const
3942 { _u.parse(evList); _v.parse(evList); }
3944 const gsFeSpace<Scalar> & rowVar()
const
3945 {
return E1::ColBlocks ? _u.rowVar() : _v.rowVar(); }
3946 const gsFeSpace<Scalar> & colVar()
const
3951 void print(std::ostream &os)
const { _u.print(os); os<<
"~"; _v.print(os); }
3955 template <
typename E1,
typename E2>
3957 collapse_expr<E1,E2> collapse( _expr<E1>
const& u, _expr<E2>
const& v)
3958 {
return collapse_expr<E1, E2>(u, v); }
3971 template <
typename E1,
typename E2,
bool = E2::ColBlocks>
3972 class frprod_expr :
public _expr<frprod_expr<E1, E2> >
3975 typedef typename E1::Scalar Scalar;
3976 enum {ScalarValued = 0, ColBlocks=E2::ColBlocks};
3977 enum { Space = (int)E1::Space + (
int)E2::Space };
3985 typename E1::Nested_t _u;
3986 typename E2::Nested_t _v;
3988 mutable gsMatrix<Scalar> res;
3992 frprod_expr(_expr<E1>
const& u, _expr<E2>
const& v)
4002 const gsMatrix<Scalar> & eval(
const index_t k)
const
4006 const index_t nb = _u.cardinality();
4007 auto A = _u.eval(k);
4008 auto B = _v.eval(k);
4010 for (
index_t i = 0; i!=nb; ++i)
4011 for (
index_t j = 0; j!=nb; ++j)
4013 (A.middleCols(i*rb,rb).array() * B.middleCols(j*rb,rb).array()).sum();
4017 index_t rows()
const {
return _u.cols() / _u.rows(); }
4018 index_t cols()
const {
return _u.cols() / _u.rows(); }
4020 void parse(gsExprHelper<Scalar> & evList)
const
4021 { _u.parse(evList); _v.parse(evList); }
4023 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4024 const gsFeSpace<Scalar> & colVar()
const {
return _v.colVar(); }
4026 void print(std::ostream &os)
const
4027 { os <<
"("; _u.print(os); os<<
" % "; _v.print(os); os<<
")";}
4037 template <
typename E1,
typename E2>
4038 class frprod_expr<E1,E2,false> :
public _expr<frprod_expr<E1, E2,false> >
4041 typedef typename E1::Scalar Scalar;
4042 enum {ScalarValued = 0, Space = E1::Space, ColBlocks= E1::ColBlocks};
4045 typename E1::Nested_t _u;
4046 typename E2::Nested_t _v;
4048 mutable gsMatrix<Scalar> res;
4052 frprod_expr(_expr<E1>
const& u, _expr<E2>
const& v)
4062 const gsMatrix<Scalar> & eval(
const index_t k)
const
4065 auto A = _u.eval(k);
4066 auto B = _v.eval(k);
4068 const index_t nb = _u.cardinality();
4070 for (
index_t i = 0; i!=nb; ++i)
4072 (A.middleCols(i*rb,rb).array() * B.array()).sum();
4076 index_t rows()
const {
return _u.cols() / _u.rows(); }
4077 index_t cols()
const {
return 1; }
4079 void parse(gsExprHelper<Scalar> & evList)
const
4080 { _u.parse(evList); _v.parse(evList); }
4083 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4084 const gsFeSpace<Scalar> & colVar()
const {
return _v.rowVar(); }
4086 void print(std::ostream &os)
const
4087 { os <<
"("; _u.print(os); os<<
" % "; _v.print(os); os<<
")";}
4093 template <
typename E1,
typename E2>
4094 class divide_expr :
public _expr<divide_expr<E1,E2> >
4096 typename E1::Nested_t _u;
4097 typename E2::Nested_t _v;
4100 typedef typename E1::Scalar Scalar;
4103 enum {ScalarValued = E1::ScalarValued, ColBlocks= E2::ColBlocks};
4104 enum {Space = E1::Space};
4106 divide_expr(_expr<E1>
const& u, _expr<E2>
const& v)
4109 GISMO_STATIC_ASSERT(E2::ScalarValued,
"The denominator needs to be scalar valued.");
4112 AutoReturn_t eval(
const index_t k)
const
4113 {
return ( _u.eval(k) / _v.eval(k) ); }
4115 index_t rows()
const {
return _u.rows(); }
4116 index_t cols()
const {
return _u.cols(); }
4118 void parse(gsExprHelper<Scalar> & evList)
const
4119 { _u.parse(evList); _v.parse(evList); }
4122 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4123 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
4125 void print(std::ostream &os)
const
4126 { os <<
"("; _u.print(os);os <<
" / ";_v.print(os);os <<
")"; }
4132 template <
typename E1>
4133 class divide_expr<E1,typename E1::Scalar>
4134 :
public _expr<divide_expr<E1,typename E1::Scalar> >
4137 typedef typename E1::Scalar Scalar;
4140 typename E1::Nested_t _u;
4144 enum {Space= E1::Space, ScalarValued = E1::ScalarValued, ColBlocks= E1::ColBlocks};
4146 divide_expr(_expr<E1>
const& u, Scalar
const c)
4149 AutoReturn_t eval(
const index_t k)
const
4150 {
return ( _u.eval(k) / _c ); }
4152 index_t rows()
const {
return _u.rows(); }
4153 index_t cols()
const {
return _u.cols(); }
4155 void parse(gsExprHelper<Scalar> & evList)
const
4156 { _u.parse(evList); }
4159 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4160 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
4162 void print(std::ostream &os)
const
4163 { os <<
"("; _u.print(os);os <<
"/"<< _c <<
")"; }
4170 template <
typename E2>
4171 class divide_expr<typename E2::Scalar,E2>
4172 :
public _expr<divide_expr<typename E2::Scalar,E2> >
4175 typedef typename E2::Scalar Scalar;
4179 typename E2::Nested_t _u;
4181 enum {Space= 0, ScalarValued = 1, ColBlocks= 0};
4183 divide_expr(Scalar
const c, _expr<E2>
const& u)
4185 { GISMO_STATIC_ASSERT(E2::ScalarValued,
"The denominator needs to be scalar valued."); }
4187 Scalar eval(
const index_t k)
const
4188 {
return ( _c / _u.val().eval(k) ); }
4190 index_t rows()
const {
return 0; }
4191 index_t cols()
const {
return 0; }
4193 void parse(gsExprHelper<Scalar> & evList)
const
4194 { _u.parse(evList); }
4197 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4198 const gsFeSpace<Scalar> & colVar()
const {
return _u.colVar(); }
4200 void print(std::ostream &os)
const
4201 { os <<
"("<< _c <<
"/";_u.print(os);os <<
")";}
4207 template <
typename E1,
typename E2>
4208 class add_expr :
public _expr<add_expr<E1, E2> >
4210 typename E1::Nested_t _u;
4211 typename E2::Nested_t _v;
4214 enum {ScalarValued = E1::ScalarValued && E2::ScalarValued,
4215 ColBlocks = E1::ColBlocks || E2::ColBlocks };
4216 enum {Space = E1::Space};
4218 typedef typename E1::Scalar Scalar;
4220 add_expr(_expr<E1>
const& u, _expr<E2>
const& v)
4224 _u.rowVar()==_v.rowVar() && _u.colVar()==_v.colVar(),
4225 "Error: adding apples and oranges (use comma instead),"
4226 " namely:\n" << _u <<
"\n"<<_v<<
4227 " \nvars:\n" << _u.rowVar().id()<<
"!="<<_v.rowVar().id() <<
", "<< _u.colVar().id()<<
"!="<<_v.colVar().id()<<
4228 " \nspaces:\n" << (int)E1::Space<<
"!="<< (
int)E2::Space
4232 mutable Temporary_t res;
4233 const Temporary_t & eval(
const index_t k)
const
4236 "Wrong dimensions "<<_u.rows()<<
"!="<<_v.rows()<<
" in + operation:\n"
4237 << _u <<
" plus \n" << _v );
4239 "Wrong dimensions "<<_u.cols()<<
"!="<<_v.cols()<<
" in + operation:\n"
4240 << _u <<
" plus \n" << _v );
4241 res = _u.eval(k) + _v.eval(k);
4245 index_t rows()
const {
return _u.rows(); }
4246 index_t cols()
const {
return _u.cols(); }
4248 void parse(gsExprHelper<Scalar> & evList)
const
4249 { _u.parse(evList); _v.parse(evList); }
4252 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
4254 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4255 const gsFeSpace<Scalar> & colVar()
const {
return _v.colVar(); }
4257 void print(std::ostream &os)
const
4258 { os <<
"("; _u.print(os);os <<
" + ";_v.print(os);os <<
")"; }
4274 template <
typename E1,
typename E2>
4275 class summ_expr :
public _expr<summ_expr<E1,E2> >
4278 typedef typename E1::Scalar Scalar;
4280 enum {Space = E1::Space, ScalarValued= 0, ColBlocks= E2::ColBlocks};
4282 summ_expr(E1
const& u, E2
const& M) : _u(u), _M(M) { }
4284 const gsMatrix<Scalar> & eval(
const index_t k)
const
4286 auto sl = _u.eval(k);
4288 auto ml = _M.eval(k);
4290 const index_t mb = _M.cardinality();
4292 GISMO_ASSERT(_M.cols()==_M.rows(),
"Matrix must be square: "<< _M.rows()<<
" x "<< _M.cols() <<
" expr: "<< _M );
4293 GISMO_ASSERT(mb==_u.cols(),
"cardinality must match vector, but card(M)="<<_M.cardinality()<<
" and cols(u)="<<_u.cols());
4295 res.setZero(mr, sr * mr);
4296 for (
index_t i = 0; i!=sr; ++i)
4297 for (
index_t t = 0; t!=mb; ++t)
4298 res.middleCols(i*mr,mr) += sl(i,t) * ml.middleCols(t*mr,mr);
4302 index_t rows()
const {
return _M.rows(); }
4303 index_t cols()
const {
return _M.rows(); }
4305 void parse(gsExprHelper<Scalar> & evList)
const
4306 { _u.parse(evList); _M.parse(evList); }
4308 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4309 const gsFeSpace<Scalar> & colVar()
const {
return gsNullExpr<Scalar>::get(); }
4311 index_t cardinality_impl()
const
4314 void print(std::ostream &os)
const
4315 { os <<
"sum("; _M.print(os); os<<
","; _u.print(os); os<<
")"; }
4318 typename E1::Nested_t _u;
4319 typename E2::Nested_t _M;
4321 mutable gsMatrix<Scalar> res;
4328 template <
typename E1,
typename E2>
4329 class sub_expr :
public _expr<sub_expr<E1, E2> >
4331 typename E1::Nested_t _u;
4332 typename E2::Nested_t _v;
4335 enum {ScalarValued = E1::ScalarValued && E2::ScalarValued,
4336 ColBlocks = E1::ColBlocks || E2::ColBlocks };
4337 enum {Space = E1::Space};
4339 typedef typename E1::Scalar Scalar;
4341 sub_expr(_expr<E1>
const& u, _expr<E2>
const& v)
4345 _u.rowVar()==_v.rowVar() && _u.colVar()==_v.colVar(),
4346 "Error: substracting apples from oranges (use comma instead),"
4347 " namely:\n" << _u <<
"\n"<<_v);
4350 mutable Temporary_t res;
4351 const Temporary_t & eval(
const index_t k)
const
4358 "Wrong dimensions "<<_u.rows()<<
"!="<<_v.rows()<<
" in - operation:\n" << _u <<
" minus \n" << _v );
4360 "Wrong dimensions "<<_u.cols()<<
"!="<<_v.cols()<<
" in - operation:\n" << _u <<
" minus \n" << _v );
4362 "Cardinality "<< _u.cardinality()<<
" != "<< _v.cardinality());
4365 res = _u.eval(k) - _v.eval(k);
4369 index_t rows()
const {
return _u.rows(); }
4370 index_t cols()
const {
return _u.cols(); }
4372 void parse(gsExprHelper<Scalar> & evList)
const
4373 { _u.parse(evList); _v.parse(evList); }
4375 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4376 const gsFeSpace<Scalar> & colVar()
const {
return _v.colVar(); }
4378 index_t cardinality_impl()
const
4381 "Cardinality "<< _u.cardinality()<<
" != "<< _v.cardinality());
4382 return _u.cardinality();
4385 void print(std::ostream &os)
const
4386 { os <<
"("; _u.print(os); os<<
" - ";_v.print(os); os <<
")";}
4392 template <
typename E>
4393 class symm_expr :
public _expr<symm_expr<E> >
4395 typename E::Nested_t _u;
4397 mutable gsMatrix<typename E::Scalar> tmp;
4399 typedef typename E::Scalar Scalar;
4401 enum { Space = (0==E::Space ? 0 : E::Space), ScalarValued= E::ScalarValued, ColBlocks= E::ColBlocks };
4403 symm_expr(_expr<E>
const& u)
4411 return tmp * tmp.transpose();
4414 index_t rows()
const {
return _u.rows(); }
4415 index_t cols()
const {
return _u.rows(); }
4417 void parse(gsExprHelper<Scalar> & evList)
const
4418 { _u.parse(evList); }
4420 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4421 const gsFeSpace<Scalar> & colVar()
const {
return _u.rowVar(); }
4423 void print(std::ostream &os)
const { os <<
"symm("; _u.print(os); os <<
")"; }
4426 template <
typename E>
4427 class symmetrize_expr :
public _expr<symmetrize_expr<E> >
4429 typename E::Nested_t _u;
4431 mutable gsMatrix<typename E::Scalar> tmp;
4433 enum { Space = (0==E::Space ? 0 : 3), ScalarValued=E::ScalarValued, ColBlocks= E::ColBlocks };
4434 typedef typename E::Scalar Scalar;
4436 symmetrize_expr(_expr<E>
const& u)
4444 return tmp + tmp.transpose();
4447 index_t rows()
const {
return _u.rows(); }
4448 index_t cols()
const {
return _u.rows(); }
4450 void parse(gsExprHelper<Scalar> & evList)
const
4451 { _u.parse(evList); }
4453 const gsFeSpace<Scalar> & rowVar()
const {
return _u.rowVar(); }
4454 const gsFeSpace<Scalar> & colVar()
const {
return _u.rowVar(); }
4455 index_t cardinality_impl()
const {
return _u.cardinality_impl(); }
4457 void print(std::ostream &os)
const { os <<
"symmetrize("; _u.print(os); os <<
")"; }
4472 EIGEN_STRONG_INLINE constMat_expr ones(
const index_t dim) {
4475 return constMat_expr(ones);
4478 EIGEN_STRONG_INLINE constMat_expr mat(
const gsMatrix<real_t> mat) {
return constMat_expr(mat); }
4485 template<
class E> EIGEN_STRONG_INLINE
4486 abs_expr<E>
abs(
const E & u) {
return abs_expr<E>(u); }
4489 template<
class E> EIGEN_STRONG_INLINE
4490 grad_expr<E>
grad(
const E & u) {
return grad_expr<E>(u); }
4493 template<
class E> EIGEN_STRONG_INLINE
4497 template<
class T> EIGEN_STRONG_INLINE
4501 template<
class T> EIGEN_STRONG_INLINE
4505 template<
class T> EIGEN_STRONG_INLINE
4508 template<
class T> EIGEN_STRONG_INLINE
4509 normal_expr<T> sn(
const gsGeometryMap<T> & u) {
return normal_expr<T>(u); }
4512 template<
class T> EIGEN_STRONG_INLINE
4515 template<
class E> EIGEN_STRONG_INLINE
4516 lapl_expr<E> lapl(
const symbol_expr<E> & u) {
return lapl_expr<E>(u); }
4518 template<
class T> EIGEN_STRONG_INLINE
4519 lapl_expr<gsFeSolution<T> > lapl(
const gsFeSolution<T> & u)
4520 {
return lapl_expr<gsFeSolution<T> >(u); }
4523 template<
class T> EIGEN_STRONG_INLINE fform2nd_expr<T>
fform2nd(
const gsGeometryMap<T> & G)
4524 {
return fform2nd_expr<T>(G); }
4527 template<
class E> EIGEN_STRONG_INLINE
4528 jac_expr<E>
jac(
const symbol_expr<E> & u) {
return jac_expr<E>(u); }
4531 template<
class T> EIGEN_STRONG_INLINE
4532 jac_expr<gsGeometryMap<T> >
jac(
const gsGeometryMap<T> & G) {
return jac_expr<gsGeometryMap<T> >(G);}
4535 template<
class T> EIGEN_STRONG_INLINE
4536 grad_expr<gsFeSolution<T> >
jac(
const gsFeSolution<T> & s) {
return grad_expr<gsFeSolution<T> >(s);}
4538 template<
class E> EIGEN_STRONG_INLINE
4539 hess_expr<E> hess(
const symbol_expr<E> & u) {
return hess_expr<E>(u); }
4542 template<
class T> EIGEN_STRONG_INLINE
4543 hess_expr<gsGeometryMap<T> > hess(
const gsGeometryMap<T> & u) {
return hess_expr<gsGeometryMap<T> >(u); }
4546 template<
class T> EIGEN_STRONG_INLINE
4547 hess_expr<gsFeSolution<T> > hess(
const gsFeSolution<T> & u) {
return hess_expr<gsFeSolution<T> >(u); }
4550 template<
class T> EIGEN_STRONG_INLINE
4551 dJacG_expr<T>
dJac(
const gsGeometryMap<T> & G) {
return dJacG_expr<T>(G); }
4554 template<
class T> EIGEN_STRONG_INLINE
4555 meas_expr<T>
meas(
const gsGeometryMap<T> & G) {
return meas_expr<T>(G); }
4558 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4559 mult_expr<E1,E2>
const operator*(_expr<E1>
const& u, _expr<E2>
const& v)
4560 {
return mult_expr<E1, E2>(u, v); }
4562 template <
typename E2> EIGEN_STRONG_INLINE
4563 mult_expr<typename E2::Scalar,E2,false>
const
4564 operator*(
typename E2::Scalar
const& u, _expr<E2>
const& v)
4565 {
return mult_expr<typename E2::Scalar, E2, false>(u, v); }
4567 template <
typename E1> EIGEN_STRONG_INLINE
4568 mult_expr<typename E1::Scalar,E1,false>
const
4569 operator*(_expr<E1>
const& v,
typename E1::Scalar
const& u)
4570 {
return mult_expr<typename E1::Scalar,E1, false>(u, v); }
4572 template <
typename E1> EIGEN_STRONG_INLINE
4573 mult_expr<typename E1::Scalar,E1,false>
const
4574 operator-(_expr<E1>
const& u)
4575 {
return mult_expr<typename E1::Scalar,E1, false>(-1, u); }
4577 template <
typename E> mult_expr<constMat_expr, E>
const
4578 operator*( gsMatrix<typename E::Scalar>
const& u, _expr<E>
const& v)
4579 {
return mult_expr<constMat_expr, E>(mat(u), v); }
4581 template <
typename E> mult_expr<E, constMat_expr>
const
4582 operator*(_expr<E>
const& u, gsMatrix<typename E::Scalar>
const& v)
4583 {
return mult_expr<E, constMat_expr>(u, mat(v) ); }
4586 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4587 frprod_expr<E1,E2>
const operator%(_expr<E1>
const& u, _expr<E2>
const& v)
4588 {
return frprod_expr<E1, E2>(u, v); }
4591 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4592 divide_expr<E1,E2>
const operator/(_expr<E1>
const& u, _expr<E2>
const& v)
4593 {
return divide_expr<E1,E2>(u, v); }
4595 template <
typename E> EIGEN_STRONG_INLINE
4596 divide_expr<E,typename E::Scalar>
const
4597 operator/(_expr<E>
const& u,
const typename E::Scalar v)
4598 {
return divide_expr<E,typename E::Scalar>(u, v); }
4600 template <
typename E> EIGEN_STRONG_INLINE
4601 divide_expr<typename E::Scalar,E>
const
4602 operator/(
const typename E::Scalar u, _expr<E>
const& v)
4603 {
return divide_expr<typename E::Scalar,E>(u, v); }
4606 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4607 add_expr<E1,E2>
const operator+(_expr<E1>
const& u, _expr<E2>
const& v)
4608 {
return add_expr<E1, E2>(u, v); }
4611 template <
typename E> EIGEN_STRONG_INLINE
4612 add_expr< E, _expr<typename E::Scalar, true> >
4614 {
return add_expr<E,_expr<typename E::Scalar>>(u, _expr<typename E::Scalar,true>(v)); }
4617 template <
typename E> EIGEN_STRONG_INLINE
4618 add_expr< E, _expr<typename E::Scalar, true> >
4620 {
return add_expr<E,_expr<typename E::Scalar>>(u, _expr<typename E::Scalar,true>(v)); }
4623 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4624 summ_expr<E1,E2>
const summ(E1
const & u, E2
const& M)
4625 {
return summ_expr<E1,E2>(u, M); }
4629 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4635 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4640 template <
typename E1,
typename E2> EIGEN_STRONG_INLINE
4641 sub_expr<E1,E2>
const operator-(_expr<E1>
const& u, _expr<E2>
const& v)
4642 {
return sub_expr<E1, E2>(u, v); }
4644 template <
typename E2> EIGEN_STRONG_INLINE
4645 sub_expr<_expr<typename E2::Scalar>,E2>
const
4646 operator-(
typename E2::Scalar
const& s, _expr<E2>
const& v)
4649 return sub_expr<_expr<typename E2::Scalar>, E2>(_expr<typename E2::Scalar>(s), v);
4655 #define GISMO_SHORTCUT_VAR_EXPRESSION(name,impl) template<class E> EIGEN_STRONG_INLINE \
4656 auto name(const E & u) -> decltype(impl) { return impl; }
4657 #define GISMO_SHORTCUT_MAP_EXPRESSION(name,impl) template<class T> EIGEN_STRONG_INLINE \
4658 auto name(const gsGeometryMap<T> & G) -> decltype(impl) { return impl; }
4659 #define GISMO_SHORTCUT_PHY_EXPRESSION(name,impl) template<class E> EIGEN_STRONG_INLINE \
4660 auto name(const E & u, const gsGeometryMap<typename E::Scalar> & G) -> decltype(impl) { return impl; }
4663 GISMO_SHORTCUT_VAR_EXPRESSION( div,
jac(u).trace() )
4664 GISMO_SHORTCUT_PHY_EXPRESSION( idiv, ijac(u,G).trace() )
4667 GISMO_SHORTCUT_MAP_EXPRESSION(unv,
nv(G).normalized() )
4669 GISMO_SHORTCUT_MAP_EXPRESSION(usn, sn(G).normalized() )
4671 GISMO_SHORTCUT_PHY_EXPRESSION(igrad, grad(u)*
jac(G).ginv() )
4672 GISMO_SHORTCUT_VAR_EXPRESSION(igrad, grad(u) )
4674 GISMO_SHORTCUT_PHY_EXPRESSION( ijac,
jac(u) *
jac(G).ginv())
4677 GISMO_SHORTCUT_PHY_EXPRESSION(ihess,
4678 jac(G).ginv().tr()*( hess(u) -
summ(igrad(u,G),hess(G)) ) *
jac(G).ginv() )
4679 GISMO_SHORTCUT_VAR_EXPRESSION(ihess, hess(u) )
4681 GISMO_SHORTCUT_PHY_EXPRESSION(ilapl, ihess(u,G).trace() )
4682 GISMO_SHORTCUT_VAR_EXPRESSION(ilapl, hess(u).trace() )
4684 GISMO_SHORTCUT_VAR_EXPRESSION(fform,
jac(u).tr()*
jac(u) )
4685 GISMO_SHORTCUT_VAR_EXPRESSION(shapeop, fform(u).inv() *
fform2nd(u) )
4687 #undef GISMO_SHORTCUT_PHY_EXPRESSION
4688 #undef GISMO_SHORTCUT_VAR_EXPRESSION
4689 #undef GISMO_SHORTCUT_MAP_EXPRESSION
EIGEN_STRONG_INLINE matrix_by_space_expr< E1, E2 > const matrix_by_space(E1 const &u, E2 const &v)
Definition: gsExpressions.h:4630
asdiag_expr< E > asDiag() const
Returns a diagonal matrix expression of the vector expression.
Definition: gsExpressions.h:309
#define MatExprType
[Include namespace]
Definition: gsThinShellUtils.h:32
mult_expr< real_t, ppart_expr< mult_expr< double, E, false > >, false > npart() const
Returns the expression's negative part.
Definition: gsExpressions.h:260
Definition: gsExpressions.h:120
Definition: gsExpressions.h:96
Definition: gsExpressions.h:137
size_t mapSize() const
Returns the total number of patch-local degrees of freedom that are being mapped. ...
Definition: gsDofMapper.h:473
Gradient of the object.
Definition: gsForwardDeclarations.h:73
EIGEN_STRONG_INLINE diag_expr< E > const diagonal(E const &u)
Get diagonal elements of matrix as a vector.
Definition: gsExpressions.h:2082
Matrix< Scalar, Dynamic, Dynamic > cramerInverse() const
Inversion for small matrices using Cramer's Rule.
Definition: gsMatrixAddons.h:55
Definition: gsExpressions.h:3044
normalized_expr< E > normalized() const
Returns the vector normalized to unit length.
Definition: gsExpressions.h:283
The functions compute Dirichlet degrees of freedom using various methods.
EIGEN_STRONG_INLINE matrix_by_space_expr_tr< E1, E2 > const matrix_by_space_tr(E1 const &u, E2 const &v)
Definition: gsExpressions.h:4636
sqNorm_expr< E > sqNorm() const
Returns the squared Euclidean norm of the expression.
Definition: gsExpressions.h:291
index_t rows() const
Returns the row-size of the expression.
Definition: gsExpressions.h:328
Definition: gsExpressions.h:2304
ppart_expr< E > ppart() const
Returns the expression's positive part.
Definition: gsExpressions.h:253
const gsFeSpace< Scalar > & rowVar() const
Definition: gsExpressions.h:358
Outward normal on the boundary.
Definition: gsForwardDeclarations.h:86
EIGEN_STRONG_INLINE tangent_expr< T > tv(const gsGeometryMap< T > &u)
The tangent boundary vector of a geometry map in 2D.
Definition: gsExpressions.h:4513
norm_expr< E > norm() const
Returns the Euclidean norm of the expression.
Definition: gsExpressions.h:279
adjugate_expr< E > adj() const
Returns the adjugate of the expression (for matrix-valued expressions)
Definition: gsExpressions.h:275
Definition: gsExpressions.h:3011
S give(S &x)
Definition: gsMemory.h:266
EIGEN_STRONG_INLINE dJacG_expr< T > dJac(const gsGeometryMap< T > &G)
The partial derivatives of the Jacobian matrix of a geometry map.
Definition: gsExpressions.h:4551
Second derivatives.
Definition: gsForwardDeclarations.h:80
#define index_t
Definition: gsConfig.h:32
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
const gsFeElement< Scalar > & _e
Reference to the element.
Definition: gsExpressions.h:856
EIGEN_STRONG_INLINE fform2nd_expr< T > fform2nd(const gsGeometryMap< T > &G)
The second fundamental form of G.
Definition: gsExpressions.h:4523
EIGEN_STRONG_INLINE frprod_expr< E1, E2 > const operator%(_expr< E1 > const &u, _expr< E2 > const &v)
Frobenious product (also known as double dot product) operator for expressions.
Definition: gsExpressions.h:4587
det_expr< E > det() const
Returns the determinant of the expression.
Definition: gsExpressions.h:287
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
EIGEN_STRONG_INLINE summ_expr< E1, E2 > const summ(E1 const &u, E2 const &M)
Matrix-summation operator for expressions.
Definition: gsExpressions.h:4624
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
Definition: gsExpressions.h:3132
Struct containing information for matrix assembly.
Definition: gsExpressions.h:63
colsum_expr< E > colSum() const
Returns the colSum of a matrix.
Definition: gsExpressions.h:321
Laplacian.
Definition: gsForwardDeclarations.h:83
const gsFeSpace< Scalar > & colVar() const
Definition: gsExpressions.h:366
void finalize()
Must be called after all boundaries and interfaces have been marked to set up the dof numbering...
Definition: gsDofMapper.cpp:240
value_expr< E > val() const
Definition: gsExpressions.h:305
temp_expr< E > temp() const
Returns an evaluation of the (sub-)expression in temporary memory.
Definition: gsExpressions.h:263
exp_expr< E > exp() const
Returns exp(expression)
Definition: gsExpressions.h:249
#define gsWarn
Definition: gsDebug.h:50
Definition: gsExpressions.h:2335
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
Normal vector of the object.
Definition: gsForwardDeclarations.h:85
inv_expr< E > const inv() const
Returns the inverse of the expression (for matrix-valued expressions)
Definition: gsExpressions.h:267
Provides declaration of Basis abstract interface.
EIGEN_STRONG_INLINE mult_expr< E1, E2 > const operator*(_expr< E1 > const &u, _expr< E2 > const &v)
Multiplication operator for expressions.
Definition: gsExpressions.h:4559
EIGEN_STRONG_INLINE onormal_expr< T > nv(const gsGeometryMap< T > &u)
The (outer pointing) boundary normal of a geometry map.
Definition: gsExpressions.h:4506
Definition: gsExpressions.h:139
Definition: gsExpressions.h:140
cb_expr< E > cb() const
Returns the puts the expression to colBlocks.
Definition: gsExpressions.h:241
tr_expr< E > tr() const
Returns the transpose of the expression.
Definition: gsExpressions.h:233
void print(std::ostream &os) const
Prints the expression as a string to os.
Definition: gsExpressions.h:194
MatExprType eval(const index_t k) const
Evaluates the expression at evaluation point indexed by k.
Definition: gsExpressions.h:229
Active function ids.
Definition: gsForwardDeclarations.h:84
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition: gsFuncData.h:23
EIGEN_STRONG_INLINE nabla_expr< T > nabla(const gsFeVariable< T > &u)
The nabla ( ) of a finite element variable.
Definition: gsExpressions.h:4502
Definition: gsExpressions.h:138
Second fundamental form.
Definition: gsForwardDeclarations.h:87
Definition: gsExpressions.h:126
The density of the measure pull back.
Definition: gsForwardDeclarations.h:76
trace_expr< E > trace() const
Returns the trace of the expression (for matrix-valued expressions)
Definition: gsExpressions.h:271
tr_expr< E, true > cwisetr() const
Returns the coordinate-wise transpose of the expression.
Definition: gsExpressions.h:237
EIGEN_STRONG_INLINE curl_expr< T > curl(const gsFeVariable< T > &u)
The curl of a finite element variable.
Definition: gsExpressions.h:4498
Definition: gsExpressions.h:2598
size_t numPatches() const
Returns the number of patches present underneath the mapper.
Definition: gsDofMapper.h:469
EIGEN_STRONG_INLINE reshape_expr< E > const reshape(E const &u, index_t n, index_t m)
Reshape an expression.
Definition: gsExpressions.h:1925
EIGEN_STRONG_INLINE dJacdc_expr< E > dJacdc(const E &u, index_t c)
The derivative of the jacobian of a geometry map with respect to a coordinate.
Definition: gsExpressions.h:4494
nabla2_expr< T > nabla2(const gsFeVariable< T > &u)
The nabla2 ( ) of a finite element variable.
Definition: gsExpressions.h:3003
Definition: gsExpressions.h:136
EIGEN_STRONG_INLINE flat_expr< E > const flat(E const &u)
Make a matrix 2x2 expression "flat".
Definition: gsExpressions.h:2031
Value of the object.
Definition: gsForwardDeclarations.h:72
EIGEN_STRONG_INLINE grad_expr< E > grad(const E &u)
The gradient of a variable.
Definition: gsExpressions.h:4490
index_t cols() const
Returns the column-size of the expression.
Definition: gsExpressions.h:332
EIGEN_STRONG_INLINE idMat_expr id(const index_t dim)
The identity matrix of dimension dim.
Definition: gsExpressions.h:4470
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
EIGEN_STRONG_INLINE meas_expr< T > meas(const gsGeometryMap< T > &G)
The measure of a geometry map.
Definition: gsExpressions.h:4555
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
EIGEN_STRONG_INLINE divide_expr< E1, E2 > const operator/(_expr< E1 > const &u, _expr< E2 > const &v)
Scalar division operator for expressions.
Definition: gsExpressions.h:4592
Gradient transformation matrix.
Definition: gsForwardDeclarations.h:77
static bool isVector()
rowSpan && !colSpan
Definition: gsExpressions.h:344
max_expr< E > max() const
Returns the rowSum of a matrix.
Definition: gsExpressions.h:313
GISMO_EXPR_VECTOR_EXPRESSION(norm, norm, 1)
Eucledian Norm.
void parse(gsExprHelper< Scalar > &evList) const
Parse the expression and discover the list of evaluation sources, also sets the required evaluation f...
Definition: gsExpressions.h:350
This object is a cache for computed values from an evaluator.
EIGEN_STRONG_INLINE add_expr< E1, E2 > const operator+(_expr< E1 > const &u, _expr< E2 > const &v)
Addition operator for expressions.
Definition: gsExpressions.h:4607
rowsum_expr< E > rowSum() const
Returns the rowSum of a matrix.
Definition: gsExpressions.h:317
Definition: gsExpressions.h:1978
EIGEN_STRONG_INLINE replicate_expr< E > const replicate(E const &u, index_t n, index_t m=1)
Replicate an expression.
Definition: gsExpressions.h:1967
Definition: gsExpressions.h:114
Defines expression precomputed_expr.
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
sign_expr< E > sgn(Scalar tolerance=0) const
Returns the sign of the expression.
Definition: gsExpressions.h:245
Definition: gsExpressions.h:3080
bool isScalar() const
Returns true iff the expression is scalar-valued.
Definition: gsExpressions.h:342
Definition: gsExpressions.h:2542
EIGEN_STRONG_INLINE jac_expr< E > jac(const symbol_expr< E > &u)
The Jacobian matrix of a FE variable.
Definition: gsExpressions.h:4528