1 /** @file gsField.h 2 3 @brief Provides declaration of the Field class. 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. Mantzaflaris 12 */ 13 14 #pragma once 15 16 #include <gsCore/gsFunctionSet.h> 17 #include <gsCore/gsGeometry.h> 18 #include <gsCore/gsMultiPatch.h> 19 #include <gsCore/gsMultiBasis.h> 20 #include <gsUtils/gsPointGrid.h> 21 22 namespace gismo 23 { 24 /** 25 * \brief A scalar of vector field defined on a m_parametric geometry. 26 * 27 * A gsField is, generally speaking, some mathematical function that is defined on a domain of interest 28 * (the name "field" is motivated by, e.g., "scalar field" or "vector field"). 29 * 30 * The gsField combines the following:\n 31 * - <b>Geometric information</b> on the domain:\n 32 * 33 * The domain can be represented as one single patch or as a 34 * collection of multiple patches (a.k.a. subdomains).\n This 35 * information is stored in a member of the gsMultiPatch class.\n 36 * 37 * - The <b>function</b> defined on the domain:\n 38 * 39 * For each patch (a.k.a. subdomain), the gsField contains a member 40 * of class gsFunction (which represents the "local field", so to 41 * say). On this, the operations of gsFunction can be carried out 42 * (e.g., function evaluation or computation of derivatives).\n 43 * Remark: The collection of patch-wise gsFunction is stored in the 44 * private member gsField::m_fields. 45 * 46 * Note that the geometry representation of a single patch can be 47 * extracted by calling the member function gsField::patch. 48 * 49 * The "local field" on a single patch can be extracted by calling gsField::function. 50 * 51 * \ingroup Core 52 */ 53 template<class T> 54 class gsField 55 { 56 57 public: 58 typedef memory::shared_ptr< gsField > Ptr;// todo: remove 59 typedef memory::unique_ptr< gsField > uPtr;// todo: remove 60 gsField()61 gsField(): m_patches(NULL) { } 62 gsField(const gsFunctionSet<T> & mp,typename gsFunctionSet<T>::Ptr fs,const bool isparam)63 gsField( const gsFunctionSet<T> & mp, 64 typename gsFunctionSet<T>::Ptr fs, 65 const bool isparam) 66 : m_patches(&mp), m_fields(fs), m_parametric(isparam) 67 { } 68 69 gsField( const gsGeometry<T> & sp, const gsFunctionSet<T> & pf, const bool isparam = false) 70 : m_patches(&sp), m_fields(memory::make_shared_not_owned(&pf)), m_parametric(isparam) 71 { } 72 gsField(const gsGeometry<T> & sp,const gsGeometry<T> & pf)73 gsField( const gsGeometry<T> & sp, const gsGeometry<T> & pf) 74 : m_patches(&sp), m_fields(memory::make_shared_not_owned(&pf)), m_parametric(true) 75 { } 76 77 gsField( const gsMultiPatch<T> & mp, const gsFunctionSet<T> & f, const bool isparam = false) 78 : m_patches(&mp), m_fields(memory::make_shared_not_owned(&f)), m_parametric(isparam) 79 { } 80 gsField(const gsMultiPatch<T> & mp,const gsMultiPatch<T> & f)81 gsField( const gsMultiPatch<T> & mp, const gsMultiPatch<T> & f) 82 : m_patches(&mp), m_fields(memory::make_shared_not_owned(&f)), m_parametric(true) 83 { } 84 85 public: 86 87 // TO DO: 88 89 // EVAL_physical_position 90 // Need to solve x = m_geometry(u) for u 91 92 // Return a point in the physical domain at parameter value u 93 /** 94 * @brief Maps points \a u from the parameter domain to the physical domain. 95 * 96 * @param[in] u Evaluation points as gsMatrix of size <em>d</em> x <em>n</em>.\n 97 * \a d denotes the dimension of the parameter domain (i.e., d = parDim()).\n 98 * \a n denotes the number of evaluation points.\n 99 * Each column of \a u corresponds to one evaluation point. 100 * @param[in] i Index of the considered patch/subdomain. 101 * @returns uPtr The <em>j</em>-th column of \a uPtr corresponds 102 * to the image of the point \a u_j (which is defined by the \a j-th column of the input parameter \a u). 103 */ 104 gsMatrix<T> point(const gsMatrix<T>& u, int i = 0) const 105 { 106 return m_patches->piece(i).eval(u); 107 } 108 109 // Return the value of the Field at parameter value u 110 // TO DO: rename to evalParam() 111 /** 112 * @brief Evaluation of the field at points \a u. 113 * 114 * @param[in] u Evaluation points as gsMatrix of size <em>d</em> x <em>n</em>.\n 115 * \a d denotes the dimension of the parameter domain (i.e., d = parDim()).\n 116 * \a n denotes the number of evaluation points.\n 117 * Each column of \a u corresponds to one evaluation point. 118 * @param[in] i Index of the considered patch/subdomain. 119 * @returns uPtr The <em>j</em>-th column of \a uPtr corresponds 120 * to the value of the field at the point \a u_j (which is defined by the \a j-th column of the input parameter \a u). 121 */ 122 gsMatrix<T> value(const gsMatrix<T>& u, int i = 0) const 123 { 124 return m_parametric 125 ? m_fields->piece(i).eval(u) 126 : m_fields->piece(i).eval( point(u, i) ); 127 } 128 129 // Return the value of the Field at physical value u 130 // TO DO: rename to evalPhys() pvalue(const gsMatrix<T> & u,int i)131 gsMatrix<T> pvalue(const gsMatrix<T>& u, int i) const 132 { 133 GISMO_ASSERT(!m_parametric, "Cannot compute physical value"); 134 return ( m_fields->piece(i).eval(u) ); 135 } 136 137 /// Computes the L2-distance between the two fields, on the physical domain 138 T distanceL2(gsField<T> const & field, int numEvals= 1000) const; 139 140 /// Computes the L2-distance between the field and a function \a func on the physical domain 141 T distanceL2(gsFunctionSet<T> const & func, 142 bool isFunc_param = false, 143 int numEvals=1000) const; 144 145 /// Computes the L2-distance between the field and a function \a 146 /// func on the physical domain, using mesh from B 147 T distanceL2(gsFunctionSet<T> const & func, 148 gsMultiBasis<T> const & B, 149 bool isFunc_param = false, 150 int numEvals=1000) const; 151 152 /// Computes the H1-seminorm of the diff. between the field and a function \a 153 /// func on the physical domain 154 T distanceH1(gsFunctionSet<T> const & func, 155 bool isFunc_param = false, 156 int = 1000) const; 157 158 /// Computes the H1-seminorm of the diff. between the field and a function \a 159 /// func on the physical domain, using mesh from B 160 T distanceH1(gsFunctionSet<T> const & func, 161 gsMultiBasis<T> const & B, 162 bool isFunc_param = false, 163 int = 1000) const; 164 165 /// Computes the H2-seminorm of the diff. between the field and a function \a 166 /// func on the physical domain, using mesh from B 167 T distanceH2(gsFunctionSet<T> const & func, 168 bool isFunc_param = false) const; 169 170 /// Computes the DG-distance between the field and a function \a 171 /// func on the physical domain 172 T distanceDG(gsFunctionSet<T> const & func, 173 bool isFunc_param = false, 174 int = 1000) const; 175 176 /// Prints the object as a string. print(std::ostream & os)177 std::ostream &print(std::ostream &os) const 178 { 179 os << ( m_parametric ? "Parameterized f" : "F") 180 << "unction field.\n Defined on " << m_patches; 181 return os; 182 } 183 184 /// \brief Returns the dimension of the parameter domain 185 /// (e.g., if the domain is a surface in three-dimensional space, it returns 2). parDim()186 short_t parDim() const { return m_patches->domainDim(); } 187 188 /// \brief Returns the dimension of the physical domain 189 /// (e.g., if the domain is a surface in three-dimensional space, it returns 3). geoDim()190 short_t geoDim() const { return m_patches->targetDim(); } 191 192 /// \brief Returns the dimension of the physical domain 193 /// (e.g., if the domain is a surface in three-dimensional space, it returns 3). dim()194 short_t dim() const { return m_fields->targetDim(); } 195 196 /// Returns the number of patches. nPatches()197 GISMO_DEPRECATED index_t nPatches() const { return m_patches->nPieces(); } 198 199 /// Returns the number of pieces. nPieces()200 int nPieces() const { return m_patches->nPieces(); } 201 geometry()202 const gsGeometry<T> & geometry() const 203 { 204 GISMO_ASSERT(dynamic_cast<const gsGeometry<T>*>(m_patches), 205 "No geometry in field. The domain is"<< *m_patches); 206 return *static_cast<const gsGeometry<T>*>(m_patches); 207 } 208 209 /// Returns gsMultiPatch containing the geometric information on the domain. patches()210 const gsMultiPatch<T> & patches() const 211 { 212 GISMO_ASSERT(dynamic_cast<const gsMultiPatch<T>*>(m_patches), 213 "No patches in field. The field domain is "<< *m_patches); 214 return *static_cast<const gsMultiPatch<T>*>(m_patches); 215 } 216 217 /// Returns the fields (defined per patch) fields()218 const gsFunctionSet<T> & fields() const { return *m_fields; } 219 220 /// Returns the gsGeometry of patch \a i. 221 const gsGeometry<T> & patch(int i=0) const 222 { 223 GISMO_ASSERT( i<nPieces(), 224 "gsField: Invalid patch index."); 225 GISMO_ASSERT(dynamic_cast<const gsGeometry<T>*>(&m_patches->piece(i)), 226 "No geometry in field. The domain is"<< m_patches->piece(i)); 227 return static_cast<const gsGeometry<T>&>(m_patches->piece(i)); 228 } 229 230 /// Returns the gsFunction of patch \a i. 231 const gsFunction<T> & function(int i=0) const 232 { 233 GISMO_ASSERT(dynamic_cast<const gsFunction<T>*>(&m_fields->piece(i)), 234 "No function in field. The domain is"<< m_patches->piece(i)); 235 return static_cast<const gsFunction<T>&>(m_fields->piece(i)); 236 } 237 238 /// Attempts to return an Isogeometric function for patch i 239 const gsGeometry<T> & igaFunction(int i=0) const 240 { 241 GISMO_ASSERT(i<m_fields->nPieces(), 242 "gsField: Invalid patch index."); 243 GISMO_ASSERT(m_parametric, 244 "Cannot get an IGA function from non-parametric field."); 245 GISMO_ASSERT(dynamic_cast<const gsGeometry<T>*>(&m_fields->piece(i)), 246 "Cannot return an igaFunction from a function of type: "<< m_fields->piece(i) ); 247 return static_cast<const gsGeometry<T> &>(m_fields->piece(i)); 248 } 249 250 /// True if the field function is defined on parametric 251 /// coordinates (i.e. the same coordinate system as the patches) isParametric()252 bool isParametric() const { return m_parametric; } 253 254 /// True if the field function is parametrized by a set of basis 255 /// functions in parametric coordinates (i.e. the same coordinate 256 /// system as the patches) isParametrized()257 bool isParametrized() const 258 { return m_parametric && dynamic_cast<const gsGeometry<T>*>(&m_fields->piece(0));} 259 260 /** \brief Returns the coefficient vector (if it exists) 261 corresponding to the function field for patch \a i. 262 263 Returns the coefficients of the field corresponding to the \a i-th 264 patch. This is only possible in the case when the field is defined 265 in terms of basis functions (ie. it derives from gsGeometry). 266 267 */ 268 const gsMatrix<T> & coefficientVector(int i=0) const 269 { 270 return igaFunction(i).coefs(); 271 } 272 273 // Data members 274 private: 275 276 /// The isogeometric field is defined on this multipatch domain 277 const gsFunctionSet<T> * m_patches; 278 279 // If there are many patches, one field per patch 280 281 /// \brief Vector containing "local fields" for each patch/subdomain. 282 /// 283 /// For each patch/subdomain, the "local field" is represented by 284 /// a gsFunction. This local field can be accessed with 285 /// gsField::function. 286 typename gsFunctionSet<T>::Ptr m_fields; 287 288 /** 289 * @brief \a True iff this is an isogeometric field. 290 * 291 * If \a m_parametric is \a true, the evaluation points for calling gsField::value have to be placed in the 292 * \a parameter domain. 293 * 294 * If \a m_parametric is \a false, then the evaluation points are in the \a physical domain. 295 * This applies to, e.g., given exact solutions which are defined on the physical domain. 296 */ 297 bool m_parametric;// True iff this is an Isogeometric field, living on parameter domain 298 299 }; // class gsField 300 301 302 /// Print (as string) operator to be used by all derived classes 303 template<class T> 304 std::ostream &operator<<(std::ostream &os, const gsField<T>& b) 305 {return b.print(os); } 306 307 308 } // namespace gismo 309 310 #ifndef GISMO_BUILD_LIB 311 #include GISMO_HPP_HEADER(gsField.hpp) 312 #endif 313