Line data Source code
1 : /*
2 : * Copyright (c) 2013-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 : * A linker that projects a given vector to the manifold of given elements.
35 : */
36 : #ifndef __H__UG__LIB_DISC__SPATIAL_DISC__PROJECT_LINKER__
37 : #define __H__UG__LIB_DISC__SPATIAL_DISC__PROJECT_LINKER__
38 :
39 : /* ug4 headers */
40 : #include "common/common.h"
41 :
42 : #include "lib_disc/reference_element/reference_mapping_provider.h"
43 : #include "linker.h"
44 :
45 : namespace ug {
46 :
47 : /**
48 : * User data projecting a given vector to the plane of degenerated elements.
49 : * The projection is computed only in the degenerated elements. For other elements,
50 : * 0 is returned.
51 : */
52 : template <int dim>
53 : class ProjectionLinker
54 : : public StdDataLinker<ProjectionLinker<dim>, MathVector<dim>, dim>
55 : {
56 : /// Base class type
57 : typedef StdDataLinker<ProjectionLinker<dim>, MathVector<dim>, dim> base_type;
58 :
59 : public:
60 :
61 : /// Constructor
62 0 : ProjectionLinker
63 : (
64 : SmartPtr<CplUserData<MathVector<dim>, dim> > spVector ///< vector to project
65 : )
66 0 : {
67 : this->set_num_input (1);
68 0 : m_spVector = spVector;
69 0 : m_spDVector = spVector.template cast_dynamic<DependentUserData<MathVector<dim>, dim> > ();
70 0 : this->set_input (0, spVector, spVector);
71 0 : }
72 :
73 : /// Constructor
74 : ProjectionLinker
75 : (
76 : MathVector<dim> vector ///< vector to project
77 : )
78 : {
79 : this->set_num_input (1);
80 : m_spVector = make_sp (new ConstUserVector<dim> (vector));
81 : m_spDVector = m_spVector.template cast_dynamic<DependentUserData<MathVector<dim>, dim> > ();
82 : this->set_input (0, m_spVector, m_spVector);
83 : }
84 :
85 : /// Returns true because without a grid function, we do not get the element to project to!
86 0 : virtual bool requires_grid_fct() const {return true;}
87 :
88 : /// Evaluation with no element is impossible
89 0 : inline void evaluate
90 : (
91 : MathVector<dim>& value,
92 : const MathVector<dim>& glob_ip,
93 : number time,
94 : int si
95 : ) const
96 : {
97 0 : UG_THROW ("ProjectionLinker: Cannot evaluate without any specification of the element!");
98 : }
99 :
100 : /// Computation only of the projection
101 : template <int refDim>
102 0 : inline void evaluate
103 : (
104 : MathVector<dim> vValue[],
105 : const MathVector<dim> vGlobIP[],
106 : number time,
107 : int si,
108 : GridObject* elem,
109 : const MathVector<dim> vCornerCoords[],
110 : const MathVector<refDim> vLocIP[],
111 : const size_t nip,
112 : LocalVector* u,
113 : const MathMatrix<refDim, dim>* vJT = NULL
114 : ) const
115 : {
116 : // 1. Get the vectors themselves
117 0 : (*m_spVector) (vValue, vGlobIP, time, si, elem, vCornerCoords, vLocIP, nip, u, vJT);
118 :
119 : // 2. Prepare the data of the element
120 0 : const ReferenceObjectID roid = elem->reference_object_id ();
121 0 : DimReferenceMapping<refDim, dim>& rMapping = ReferenceMappingProvider::get<refDim, dim> (roid);
122 0 : rMapping.update (vCornerCoords);
123 :
124 : // 3. At every integration point
125 0 : for (size_t i = 0; i < nip; i++)
126 : {
127 : // 3a. Get the Jacobian
128 : MathMatrix<dim, refDim> J;
129 0 : rMapping.jacobian (J, vLocIP[i]);
130 : // 3b. Project the vector to the subspace spanned by the columns of the Jacobian
131 0 : OrthogProjectVec (vValue[i], J);
132 : }
133 0 : }
134 :
135 : /// Computation of the projection and its derivatives
136 : template <int refDim>
137 0 : void eval_and_deriv
138 : (
139 : MathVector<dim> vValue[],
140 : const MathVector<dim> vGlobIP[],
141 : number time,
142 : int si,
143 : GridObject* elem,
144 : const MathVector<dim> vCornerCoords[],
145 : const MathVector<refDim> vLocIP[],
146 : const size_t nip,
147 : LocalVector* u,
148 : bool bDeriv,
149 : int s,
150 : std::vector<std::vector<MathVector<dim> > > vvvDeriv[],
151 : const MathMatrix<refDim, dim>* vJT = NULL
152 : ) const
153 : {
154 : // 1. Get the vectors themselves
155 0 : const MathVector<dim>* vVector = m_spVector->values (s);
156 :
157 : // 2a. Prepare the data of the element
158 0 : const ReferenceObjectID roid = elem->reference_object_id ();
159 0 : DimReferenceMapping<refDim, dim>& rMapping = ReferenceMappingProvider::get<refDim, dim> (roid);
160 0 : rMapping.update (vCornerCoords);
161 :
162 : // 2b. Check if we should compute the derivatives
163 0 : if (this->zero_derivative ())
164 : bDeriv = false;
165 : else
166 0 : this->set_zero (vvvDeriv, nip);
167 :
168 : // 3. At every integration point
169 0 : for (size_t i = 0; i < nip; i++)
170 : {
171 : // 3a. Get the Jacobian
172 : MathMatrix<dim, refDim> J;
173 0 : rMapping.jacobian (J, vLocIP[i]);
174 : // 3b. Project the vector to the subspace spanned by the columns of the Jacobian
175 0 : vValue[i] = vVector[i];
176 0 : OrthogProjectVec (vValue[i], J);
177 : // 3c. Project the derivatives in the same way
178 0 : if (! bDeriv) continue;
179 0 : for (size_t fct = 0; fct < m_spDVector->num_fct(); fct++)
180 : {
181 : const MathVector<dim>* vDVector = m_spDVector->deriv (s, i, fct);
182 : const size_t c_fct = this->input_common_fct (0, fct);
183 0 : for (size_t sh = 0; sh < this->num_sh (c_fct); sh++)
184 : {
185 0 : vvvDeriv[i][c_fct][sh] = vDVector [sh];
186 0 : OrthogProjectVec (vvvDeriv[i][c_fct][sh], J);
187 : }
188 : }
189 : }
190 0 : }
191 :
192 : private:
193 :
194 : /// data to project
195 : SmartPtr<CplUserData<MathVector<dim>, dim> > m_spVector;
196 : SmartPtr<DependentUserData<MathVector<dim>, dim> > m_spDVector;
197 : };
198 :
199 : } // end namespace ug
200 :
201 : #endif // __H__UG__LIB_DISC__SPATIAL_DISC__PROJECT_LINKER__
202 :
203 : /* End of File */
|
