Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
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 : #include "fe_geom.h"
34 : #include "lib_disc/domain_traits.h"
35 : #include "lib_disc/common/geometry_util.h"
36 : #include "lib_disc/quadrature/quadrature_provider.h"
37 :
38 : namespace ug{
39 :
40 : ////////////////////////////////////////////////////////////////////////////////
41 : // DimFEGeometry
42 : ////////////////////////////////////////////////////////////////////////////////
43 :
44 : template <int TWorldDim, int TRefDim>
45 0 : DimFEGeometry<TWorldDim,TRefDim>::
46 : DimFEGeometry() :
47 0 : m_roid(ROID_UNKNOWN), m_quadOrder(0),
48 : m_lfeID(),
49 0 : m_vIPLocal(NULL), m_vQuadWeight(NULL)
50 0 : {}
51 :
52 : template <int TWorldDim, int TRefDim>
53 0 : DimFEGeometry<TWorldDim,TRefDim>::
54 : DimFEGeometry(size_t order, LFEID lfeid) :
55 0 : m_roid(ROID_UNKNOWN), m_quadOrder(order), m_lfeID(lfeid),
56 0 : m_vIPLocal(NULL), m_vQuadWeight(NULL)
57 0 : {}
58 :
59 : template <int TWorldDim, int TRefDim>
60 0 : DimFEGeometry<TWorldDim,TRefDim>::
61 : DimFEGeometry(ReferenceObjectID roid, size_t order, LFEID lfeid) :
62 0 : m_roid(roid), m_quadOrder(order), m_lfeID(lfeid),
63 0 : m_vIPLocal(NULL), m_vQuadWeight(NULL)
64 0 : {}
65 :
66 : template <int TWorldDim, int TRefDim>
67 : void
68 0 : DimFEGeometry<TWorldDim,TRefDim>::
69 : update_local(ReferenceObjectID roid, const LFEID& lfeID, size_t orderQuad)
70 : {
71 : // remember current setting
72 0 : m_roid = roid;
73 0 : m_lfeID = lfeID;
74 0 : m_quadOrder = orderQuad;
75 :
76 : // request for quadrature rule
77 : try{
78 : const QuadratureRule<dim>& quadRule
79 0 : = QuadratureRuleProvider<dim>::get(roid, orderQuad);
80 :
81 : // copy quad informations
82 0 : m_nip = quadRule.size();
83 0 : m_vIPLocal = quadRule.points();
84 0 : m_vQuadWeight = quadRule.weights();
85 :
86 0 : }UG_CATCH_THROW("FEGeometry::update: Quadrature Rule error.");
87 :
88 : // resize for number of integration points
89 0 : m_vIPGlobal.resize(m_nip);
90 0 : m_vJTInv.resize(m_nip);
91 0 : m_vDetJ.resize(m_nip);
92 :
93 0 : m_vvGradGlobal.resize(m_nip);
94 0 : m_vvGradLocal.resize(m_nip);
95 0 : m_vvShape.resize(m_nip);
96 :
97 : // request for trial space
98 : try{
99 : const LocalShapeFunctionSet<dim>& lsfs
100 0 : = LocalFiniteElementProvider::get<dim>(roid, m_lfeID);
101 :
102 : // copy shape infos
103 0 : m_nsh = lsfs.num_sh();
104 :
105 : // resize for number of shape functions
106 0 : for(size_t ip = 0; ip < m_nip; ++ip)
107 : {
108 0 : m_vvGradGlobal[ip].resize(m_nsh);
109 0 : m_vvGradLocal[ip].resize(m_nsh);
110 0 : m_vvShape[ip].resize(m_nsh);
111 : }
112 :
113 : // get all shapes by on call
114 0 : for(size_t ip = 0; ip < m_nip; ++ip)
115 : {
116 0 : lsfs.shapes(&(m_vvShape[ip][0]), m_vIPLocal[ip]);
117 0 : lsfs.grads(&(m_vvGradLocal[ip][0]), m_vIPLocal[ip]);
118 : }
119 :
120 0 : }UG_CATCH_THROW("FEGeometry::update: Shape Function error.");
121 0 : }
122 :
123 : template <int TWorldDim, int TRefDim>
124 : void
125 0 : DimFEGeometry<TWorldDim,TRefDim>::
126 : update(GridObject* pElem, const MathVector<worldDim>* vCorner,
127 : const LFEID& lfeID, size_t orderQuad)
128 : {
129 : // check if same element
130 0 : if(pElem == m_pElem) return;
131 0 : else m_pElem = pElem;
132 :
133 : // get reference element type
134 0 : ReferenceObjectID roid = pElem->reference_object_id();
135 :
136 : // if already prepared for this roid, skip update of local values
137 0 : if(roid != m_roid || lfeID != m_lfeID || (int)orderQuad != m_quadOrder)
138 0 : update_local(roid, lfeID, orderQuad);
139 :
140 : // get reference element mapping
141 : try{
142 : DimReferenceMapping<dim, worldDim>& map
143 : = ReferenceMappingProvider::get<dim, worldDim>(roid, vCorner);
144 :
145 : // compute global integration points
146 0 : map.local_to_global(&(m_vIPGlobal[0]), &(m_vIPLocal[0]), m_nip);
147 :
148 : // compute transformation inverse and determinate at ip
149 0 : map.jacobian_transposed_inverse(&(m_vJTInv[0]), &(m_vDetJ[0]),
150 0 : &(m_vIPLocal[0]), m_nip);
151 :
152 : // compute global gradients
153 0 : for(size_t ip = 0; ip < m_nip; ++ip)
154 0 : for(size_t sh = 0; sh < m_nsh; ++sh)
155 : MatVecMult(m_vvGradGlobal[ip][sh],
156 : m_vJTInv[ip], m_vvGradLocal[ip][sh]);
157 :
158 0 : }UG_CATCH_THROW("FEGeometry::update: Reference Mapping error.");
159 : }
160 :
161 : template <int TWorldDim, int TRefDim>
162 : void
163 0 : DimFEGeometry<TWorldDim,TRefDim>::
164 : update_boundary_faces(GridObject* pElem,
165 : const MathVector<worldDim>* vCornerCoords,
166 : size_t quadOrder,
167 : const ISubsetHandler* ish)
168 : {
169 : typedef typename domain_traits<dim>::side_type Side;
170 :
171 : // get reference element type
172 0 : const ReferenceObjectID roid = pElem->reference_object_id();
173 :
174 : // get grid
175 0 : Grid& grid = *(ish->grid());
176 :
177 : // vector of subset indices of side
178 : std::vector<int> vSubsetIndex;
179 :
180 : // get subset indices for sides (i.e. edge in 2d, faces in 3d)
181 : std::vector<Side*> vSide;
182 0 : CollectAssociated(vSide, grid, pElem);
183 0 : vSubsetIndex.resize(vSide.size());
184 0 : for(size_t i = 0; i < vSide.size(); ++i)
185 0 : vSubsetIndex[i] = ish->get_subset_index(vSide[i]);
186 :
187 : // get reference element mapping
188 : try{
189 : DimReferenceMapping<dim, worldDim>& rMapping
190 0 : = ReferenceMappingProvider::get<dim, worldDim>(roid);
191 :
192 : // update reference mapping
193 0 : rMapping.update(vCornerCoords);
194 :
195 : const DimReferenceElement<dim>& rRefElem
196 : = ReferenceElementProvider::get<dim>(roid);
197 :
198 : const LocalShapeFunctionSet<dim>& rTrialSpace =
199 0 : LocalFiniteElementProvider::get<dim>(m_roid, m_lfeID);
200 :
201 : // loop requested subset
202 : typename std::map<int, std::vector<BF> >::iterator it;
203 0 : for (it=m_mapVectorBF.begin() ; it != m_mapVectorBF.end(); ++it)
204 : {
205 : // get subset index
206 0 : const int bndIndex = (*it).first;
207 :
208 : // get vector of BF for element
209 0 : std::vector<BF>& vBF = (*it).second;
210 :
211 : // clear vector
212 : vBF.clear();
213 :
214 : // loop sides of element
215 0 : for(size_t side = 0; side < vSubsetIndex.size(); ++side)
216 : {
217 : // skip non boundary sides
218 0 : if(vSubsetIndex[side] != bndIndex) continue;
219 :
220 0 : vBF.resize(vBF.size()+1);
221 : BF& bf = vBF.back();
222 :
223 : const ReferenceObjectID sideRoid = rRefElem.roid(dim-1,side);
224 :
225 0 : std::vector<MathVector<worldDim> > vSideCorner(rRefElem.num(dim-1, side, 0));
226 0 : std::vector<MathVector<dim> > vLocalSideCorner(rRefElem.num(dim-1, side, 0));
227 0 : for(size_t co = 0; co < vSideCorner.size(); ++co){
228 0 : vSideCorner[co] = vCornerCoords[rRefElem.id(dim-1, side, 0, co)];
229 : vLocalSideCorner[co] = rRefElem.corner(rRefElem.id(dim-1, side, 0, co));
230 : }
231 :
232 : const QuadratureRule<dim-1>& rSideQuadRule
233 0 : = QuadratureRuleProvider<dim-1>::get(sideRoid, quadOrder);
234 :
235 : // normal on scvf
236 0 : ElementNormal<worldDim>(sideRoid, bf.Normal, &vSideCorner[0]);
237 :
238 : // compute volume
239 0 : bf.Vol = VecTwoNorm(bf.Normal);
240 :
241 : // compute local integration points
242 0 : bf.vWeight = rSideQuadRule.weights();
243 0 : bf.nip = rSideQuadRule.size();
244 0 : bf.vLocalIP.resize(bf.nip);
245 0 : bf.vGlobalIP.resize(bf.nip);
246 :
247 : DimReferenceMapping<dim-1, dim>& map
248 : = ReferenceMappingProvider::get<dim-1, dim>(sideRoid, vLocalSideCorner);
249 0 : for(size_t ip = 0; ip < rSideQuadRule.size(); ++ip)
250 0 : map.local_to_global(bf.vLocalIP[ip], rSideQuadRule.point(ip));
251 :
252 : // compute global integration points
253 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
254 0 : rMapping.local_to_global(bf.vGlobalIP[ip], bf.vLocalIP[ip]);
255 :
256 0 : bf.vJtInv.resize(bf.nip);
257 0 : bf.vDetJ.resize(bf.nip);
258 :
259 0 : rMapping.jacobian_transposed_inverse(&bf.vJtInv[0], &bf.vLocalIP[0], bf.num_ip());
260 : DimReferenceMapping<dim-1, worldDim>& map2
261 : = ReferenceMappingProvider::get<dim-1, worldDim>(sideRoid, vSideCorner);
262 0 : map2.sqrt_gram_det(&bf.vDetJ[0], rSideQuadRule.points(), bf.num_ip());
263 :
264 0 : bf.nsh = rTrialSpace.num_sh();
265 0 : bf.vvShape.resize(bf.nip);
266 0 : bf.vvLocalGrad.resize(bf.nip);
267 0 : bf.vvGlobalGrad.resize(bf.nip);
268 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
269 : {
270 0 : bf.vvShape[ip].resize(bf.nsh);
271 0 : bf.vvLocalGrad[ip].resize(bf.nsh);
272 0 : bf.vvGlobalGrad[ip].resize(bf.nsh);
273 : }
274 :
275 : // compute shapes and gradients
276 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
277 : {
278 0 : rTrialSpace.shapes(&(bf.vvShape[ip][0]), bf.local_ip(ip));
279 0 : rTrialSpace.grads(&(bf.vvLocalGrad[ip][0]), bf.local_ip(ip));
280 : }
281 :
282 : // compute global gradient
283 0 : for(size_t ip = 0; ip < bf.num_ip(); ++ip)
284 0 : for(size_t sh = 0 ; sh < bf.num_sh(); ++sh)
285 : MatVecMult(bf.vvGlobalGrad[ip][sh],
286 : bf.vJtInv[ip], bf.vvLocalGrad[ip][sh]);
287 :
288 : }
289 : }
290 : }
291 0 : UG_CATCH_THROW("DimFEGeometry: update failed.");
292 :
293 0 : }
294 : ////////////////////////////////////////////////////////////////////////////////
295 : // explicit instantiations
296 : ////////////////////////////////////////////////////////////////////////////////
297 :
298 : template class DimFEGeometry<1, 1>;
299 : template class DimFEGeometry<2, 1>;
300 : template class DimFEGeometry<3, 1>;
301 :
302 : template class DimFEGeometry<2, 2>;
303 : template class DimFEGeometry<3, 2>;
304 :
305 : template class DimFEGeometry<3, 3>;
306 :
307 : } // end namespace ug
|
