1 /** @file gsFunctionWithDerivatives.h
2 
3     @brief A function with explicitly set derivatives
4 
5     This file is part of the G+Smo library.
6 
7     This Source Code Form is subject to the terms of the Mozilla Public
8     License, v. 2.0. If a copy of the MPL was not distributed with this
9     file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11     Author(s): A. Bressan
12 **/
13 
14 #include <gsCore/gsFunction.h>
15 
16 
17 #pragma once
18 
19 namespace gismo {
20 
21 template <typename T=real_t>
22 class gsFunctionWithDerivatives : public gsFunction<T>
23 {
24 protected:
25     const gsFunction<T> *m_values;
26     const gsFunction<T> *m_derivs;
27     const gsFunction<T> *m_deriv2;
28 public:
29     /// Shared pointer for gsFunctionWithDerivatives
30     typedef memory::shared_ptr< gsFunctionWithDerivatives > Ptr;
31 
32     /// Unique pointer for gsFunctionWithDerivatives
33     typedef memory::unique_ptr< gsFunctionWithDerivatives > uPtr;
34 
gsFunctionWithDerivatives()35     gsFunctionWithDerivatives()
36     : m_values(NULL),
37       m_derivs(NULL),
38       m_deriv2(NULL)
39     { }
40 
gsFunctionWithDerivatives(const gsFunction<T> & func,const gsFunction<T> & deriv)41     gsFunctionWithDerivatives(const gsFunction<T> &func, const gsFunction<T> &deriv)
42     : m_values(&func),
43       m_derivs(&deriv),
44       m_deriv2(NULL)
45     {
46         GISMO_ASSERT(checkDimensions(), "Dimensions do not fit");
47     }
48 
gsFunctionWithDerivatives(const gsFunction<T> & func,const gsFunction<T> & deriv,const gsFunction<T> & deriv2)49     gsFunctionWithDerivatives(const gsFunction<T> &func,
50                               const gsFunction<T> &deriv,
51                               const gsFunction<T> &deriv2)
52     : m_values(&func  ),
53       m_derivs(&deriv ),
54       m_deriv2(&deriv2)
55     {
56         GISMO_ASSERT(checkDimensions(), "Dimensions do not fit");
57     }
58 
eval_into(const gsMatrix<T> & u,gsMatrix<T> & result)59     void eval_into(const gsMatrix<T> &u, gsMatrix<T> &result) const
60     {
61         GISMO_ASSERT(m_derivs,"Not initialized");
62         m_values->eval_into(u,result);
63     }
64 
deriv_into(const gsMatrix<T> & u,gsMatrix<T> & result)65     void deriv_into(const gsMatrix<T> &u, gsMatrix<T> &result) const
66     {
67         GISMO_ASSERT(m_derivs,"No first derivative available");
68         m_derivs->eval_into(u,result);
69     }
70 
deriv2_into(const gsMatrix<T> & u,gsMatrix<T> & result)71     void deriv2_into(const gsMatrix<T> &u, gsMatrix<T> &result) const
72     {
73         GISMO_ASSERT(m_deriv2,"No second derivative available");
74         m_deriv2->eval_into(u,result);
75     }
76 
laplacian(const gsMatrix<T> & u)77     gsMatrix<T> laplacian(const gsMatrix<T> &u) const
78     {
79         gsMatrix<T>  secDer;
80         deriv2_into(u,secDer);
81         return secDer.topRows(m_values->domainDim()).colwise().sum();
82     }
83 
getDeriv()84     const gsFunction<T> & getDeriv () const
85     {
86         return *m_derivs;
87     }
getDeriv2()88     const gsFunction<T> & getDeriv2 () const
89     {
90         return *m_deriv2;
91     }
92 
piece(const index_t)93     const gsFunctionWithDerivatives & piece(const index_t) const
94     {
95         // same on all pieces
96         return *this;
97     }
98 
GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(gsFunction<T>,clone)99     GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(gsFunction<T>, clone)
100 
101     short_t targetDim () const
102     {
103         return m_values->targetDim();
104     }
domainDim()105     short_t domainDim () const
106     {
107         return m_values->domainDim();
108     }
109 
110 
111     /*
112     bool check()
113     {
114         if (m_derivs)
115         {
116 
117         }
118         if (m_deriv2)
119         {
120 
121         }
122     }
123     */
124 
125 private:
126 
127     gsFunctionWithDerivatives(const gsFunctionWithDerivatives<T> &func);
128 
checkDimensions()129     bool checkDimensions()
130     {
131         bool ok=true;
132 
133         const short_t parDim=m_values->domainDim();
134         const short_t tarDim=m_values->targetDim();
135 
136         if (m_derivs)
137         {
138             ok = ok &&  m_derivs->domainDim() == parDim;
139             ok = ok &&  m_derivs->targetDim() == parDim*tarDim;
140         }
141         if (m_deriv2)
142         {
143             ok = ok &&  m_deriv2->domainDim() == parDim;
144             ok = ok &&  m_deriv2->targetDim() == (parDim+1)*parDim*tarDim/2;
145         }
146         return ok;
147     }
148 };
149 
150 
151 }
152