Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Dmitry Logashenko
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : /*
34 : * Classes for user data that compute normals or project vectors to low-dimensional
35 : * subsets. Note that these classes are mainly implemented for visualization purposes
36 : * and cannot be used for coupling.
37 : */
38 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__ONC_USER_DATA__
39 : #define __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__ONC_USER_DATA__
40 :
41 : #include <vector>
42 :
43 : // ug4 headers
44 : #include "common/common.h"
45 : #include "common/math/ugmath.h"
46 : #include "lib_grid/grid_objects/grid_dim_traits.h"
47 : #include "lib_disc/domain_traits.h"
48 : #include "lib_disc/domain_util.h"
49 : #include "lib_disc/reference_element/reference_mapping_provider.h"
50 : #include "lib_disc/spatial_disc/user_data/std_user_data.h"
51 :
52 : namespace ug {
53 :
54 : /// Projection to the outer normal for a given vector user data
55 : /**
56 : * This class implements the UserData for evaluation of the normal component of vectors
57 : * computed by the other user data object. The normal component is evaluation only for
58 : * WDim-1-dimensional subsets - sides of WDim-dimensional elements.
59 : */
60 : template <typename TDomain>
61 : class OutNormCmp
62 : : public StdUserData<OutNormCmp<TDomain>, MathVector<TDomain::dim>, TDomain::dim, void, UserData<MathVector<TDomain::dim>, TDomain::dim, void> >
63 : {
64 : /// the world dimension
65 : static const int dim = TDomain::dim;
66 :
67 : /// the domain type
68 : typedef TDomain domain_type;
69 :
70 : /// the grid type
71 : typedef typename TDomain::grid_type grid_type;
72 :
73 : /// the original data type
74 : typedef MathVector<dim> vec_type;
75 :
76 : /// the original data
77 : SmartPtr<UserData<vec_type, dim, void> > m_spData;
78 :
79 : /// the domain
80 : SmartPtr<domain_type> m_spDomain;
81 :
82 : /// the subset group for the inner part
83 : SubsetGroup m_ssGrp;
84 :
85 : public:
86 :
87 : /// constructor
88 0 : OutNormCmp
89 : (
90 : SmartPtr<domain_type> spDomain, ///< the domain
91 : SmartPtr<UserData<vec_type, dim, void> > spData, ///< the original data
92 : const char * ss_names ///< the subsets
93 : )
94 0 : : m_spData (spData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
95 : {
96 : // Parse the subset names:
97 : std::vector<std::string> vssNames;
98 : try
99 : {
100 0 : TokenizeString (ss_names, vssNames);
101 0 : for (size_t k = 0; k < vssNames.size (); k++)
102 0 : RemoveWhitespaceFromString (vssNames [k]);
103 : m_ssGrp.clear ();
104 0 : m_ssGrp.add (vssNames);
105 0 : } UG_CATCH_THROW ("SubsetIndicatorUserData: Failed to parse subset names.");
106 0 : }
107 :
108 : /// constructor
109 0 : OutNormCmp
110 : (
111 : SmartPtr<domain_type> spDomain, ///< the domain
112 : SmartPtr<UserData<vec_type, dim, void> > spData ///< the original data
113 : )
114 0 : : m_spData (spData), m_spDomain (spDomain), m_ssGrp (spDomain->subset_handler ())
115 0 : {}
116 :
117 : /// Indicator functions are discontinuous
118 0 : virtual bool continuous () const {return false;}
119 :
120 : /// Returns true to get the grid element in the evaluation routine
121 0 : virtual bool requires_grid_fct () const {return m_spData->requires_grid_fct ();}
122 : // Remark: Note that actually we have not enought data for that. This dependence is neglected.
123 :
124 : /// Evaluator
125 : template <int refDim>
126 0 : inline void evaluate
127 : (
128 : vec_type vValue[],
129 : const MathVector<dim> vGlobIP [],
130 : number time,
131 : int si,
132 : GridObject * elem,
133 : const MathVector<dim> vCornerCoords [],
134 : const MathVector<refDim> vLocIP [],
135 : const size_t nip,
136 : LocalVector * u,
137 : const MathMatrix<refDim, dim> * vJT = NULL
138 : ) const
139 : {
140 : if (refDim == dim)
141 : {
142 0 : (* m_spData) (vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
143 0 : return;
144 : }
145 : if (refDim != dim - 1)
146 : {
147 0 : UG_THROW ("OutNormCmp: Wrong dimensionality of the subset.");
148 : }
149 :
150 : // Get the full-dim. element associated with the given side
151 : typedef typename grid_dim_traits<dim>::grid_base_object fd_elem_type;
152 : fd_elem_type * fd_elem = NULL;
153 : int fd_si = -1;
154 :
155 : typedef typename Grid::traits<fd_elem_type>::secure_container fd_elem_secure_container_type;
156 : fd_elem_secure_container_type fd_elem_list;
157 0 : grid_type& grid = * (grid_type*) (m_spDomain->grid().get ());
158 0 : grid.associated_elements (fd_elem_list, elem);
159 :
160 0 : if (m_ssGrp.empty ()) // no subsets specified; we assume, elem should be the bnd of the whole domain
161 : {
162 0 : if (fd_elem_list.size () == 1)
163 : {
164 : fd_elem = fd_elem_list[0];
165 0 : fd_si = m_spDomain->subset_handler()->get_subset_index (fd_elem);
166 : }
167 : // otherwise this is no boundary of the domain and we return 0
168 : }
169 : else
170 : {
171 0 : for (size_t k = 0; k < fd_elem_list.size (); k++)
172 : {
173 : fd_elem_type * e = fd_elem_list[k];
174 0 : int e_si = m_ssGrp.subset_handler()->get_subset_index (e);
175 0 : if (m_ssGrp.contains (e_si))
176 : {
177 0 : if (fd_elem != NULL) // i.e. e is already the second one
178 : { // this is no boundary of the subset; return 0
179 : fd_elem = NULL;
180 : break;
181 : }
182 : fd_elem = e;
183 : fd_si = e_si;
184 : }
185 : }
186 : }
187 0 : if (fd_elem == NULL)
188 : { // this is no bondary, return 0
189 0 : for (size_t ip = 0; ip < nip; ip++) vValue[ip] = 0;
190 : return;
191 : }
192 :
193 : // Get the coordinates of the corners of the full-dim. element
194 0 : std::vector<MathVector<dim> > fd_elem_corner_coords (domain_traits<dim>::MaxNumVerticesOfElem);
195 0 : CollectCornerCoordinates (fd_elem->base_object_id (), fd_elem_corner_coords,
196 : *fd_elem, *m_spDomain, true);
197 :
198 : // Convert the local ip coordinates (use the global coordinates as they are the same)
199 0 : const ReferenceObjectID fd_roid = fd_elem->reference_object_id ();
200 : DimReferenceMapping<dim, dim> & fd_map
201 : = ReferenceMappingProvider::get<dim, dim> (fd_roid, fd_elem_corner_coords);
202 0 : std::vector<MathVector<dim> > fd_loc_ip (nip);
203 0 : for (size_t ip = 0; ip < nip; ip++)
204 0 : fd_map.global_to_local(fd_loc_ip[ip], vGlobIP[ip]);
205 :
206 : // Call the original UserData, get the vectors
207 0 : (* m_spData) (&(vValue[0]), vGlobIP, time, fd_si, fd_elem, &(fd_elem_corner_coords[0]),
208 : &(fd_loc_ip[0]), nip, u);
209 : //TODO: Note that u is here a dummy parameter: We do not have enough data for it.
210 :
211 : // Project the vectors
212 0 : const ReferenceObjectID roid = elem->reference_object_id ();
213 : DimReferenceMapping<refDim, dim> & map
214 : = ReferenceMappingProvider::get<refDim, dim> (roid, vCornerCoords);
215 0 : for (size_t ip = 0; ip < nip; ip++)
216 : {
217 : MathMatrix<dim, refDim> J;
218 0 : MathVector<dim> p (vValue[ip]);
219 :
220 0 : map.jacobian (J, vLocIP[ip]);
221 0 : OrthogProjectVec (p, J);
222 : vValue[ip] -= p;
223 : }
224 : };
225 :
226 : /// This function should not be used
227 0 : void operator()
228 : (
229 : vec_type & vValue,
230 : const MathVector<dim> & globIP,
231 : number time,
232 : int si
233 : )
234 : const
235 : {
236 0 : UG_THROW("OutNormCmp: Element required for evaluation, but not passed. Cannot evaluate.");
237 : }
238 :
239 : /// This function should not be used
240 0 : void operator()
241 : (
242 : vec_type vValue [],
243 : const MathVector<dim> vGlobIP [],
244 : number time,
245 : int si,
246 : const size_t nip
247 : ) const
248 : {
249 0 : UG_THROW("OutNormCmp: Element required for evaluation, but not passed. Cannot evaluate.");
250 : }
251 :
252 : };
253 :
254 : } // end namespace ug
255 :
256 : #endif // __H__UG__LIB_DISC__SPATIAL_DISC__USER_DATA__ONC_USER_DATA__
257 :
258 : /* End of File */
|
