Line data Source code
1 : /*
2 : * Copyright (c) 2015: G-CSC, Goethe University Frankfurt
3 : * Author: Michael Lampe
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 : * Singular sources and sinks for the FV discretizations: Implementation of the functions
35 : */
36 :
37 : // ug4 headers
38 : #include "lib_grid/algorithms/ray_element_intersection_util.h"
39 : #include "lib_disc/reference_element/reference_element_util.h"
40 :
41 : namespace ug{
42 :
43 : /**
44 : * Test whether a source/sink point corresponds to a given corner of the element.
45 : * Note that \c TElem must specify the dimensionality of the element. Thus,
46 : * \c GridObject is not a proper type for \c TElem. It should be at least
47 : * \c Edge, \c Face, \c Volume or derived classes.
48 : */
49 : template <int dim, typename TData>
50 : template <typename TElem, typename TAAPos, typename TFVGeom>
51 0 : bool FVPointSourceOrSink<dim, TData>::corresponds_to
52 : (
53 : TElem* elem, ///< the element
54 : Grid& grid, ///< the grid
55 : TAAPos& aaPos, ///< position of the vertices
56 : const TFVGeom& geo, ///< FV geometry (initialized for 'elem')
57 : size_t co ///< corner to get the contribution for
58 : )
59 : {
60 : // restrict point to element
61 0 : if (!ContainsPoint(elem, point, aaPos))
62 : return false;
63 :
64 : // restrict point to scv
65 0 : for (size_t i = 0; i < geo.num_scvf(); i++)
66 : {
67 : const typename TFVGeom::SCVF& scvf = geo.scvf(i);
68 : MathVector<dim> d;
69 : VecSubtract(d, point, scvf.global_ip());
70 0 : if (scvf.from() == co)
71 : {
72 0 : if (VecDot(d, scvf.normal()) >= 0.0)
73 0 : return false;
74 : }
75 0 : else if (scvf.to() == co)
76 : {
77 0 : if (VecDot(d, scvf.normal()) < 0.0)
78 : return false;
79 : }
80 : }
81 : return true;
82 : };
83 :
84 : /**
85 : * Test whether a line source/sink corresponds to a given corner of an element.
86 : * If yes, the function returns true and sets ls and le to the beginning and the
87 : * end of the subsegment of the source/sink.
88 : *
89 : * Note that \c TElem must specify the dimensionality of the element. Thus,
90 : * \c GridObject is not a proper type for \c TElem. It should be at least
91 : * \c Edge, \c Face, \c Volume or derived classes.
92 : *
93 : * \todo This function works only for there is only one scv per corner. This is not true for pyramids!
94 : */
95 : template <int dim, typename TData>
96 : template <typename TElem, typename TAAPos, typename TFVGeom>
97 0 : bool FVLineSourceOrSink<dim, TData>::corresponds_to
98 : (
99 : TElem* elem, ///< [in] the element
100 : Grid& grid, ///< [in] the grid
101 : TAAPos& aaPos, ///< [in] position of the vertices
102 : const TFVGeom& geo, ///< [in] FV geometry (initialized for 'elem')
103 : size_t co, ///< [in] corner to get the contribution for
104 : MathVector<dim>& ls, ///< [out] beginning of the subsegment
105 : MathVector<dim>& le ///< [out] end of the subsegment
106 : )
107 : {
108 : // get the dimensionality of the element
109 0 : int elem_dim = ReferenceElementDimension(elem->reference_object_id());
110 :
111 : // restrict line segment to element (so that ls and le are the intersection points with its sides)
112 : ls = point1;
113 : le = point2;
114 : MathVector<dim> dir;
115 : number lambda_min, lambda_max;
116 : VecSubtract(dir, le, ls);
117 0 : if (!RayElementIntersection(lambda_min, lambda_max, ls, dir, elem, grid, aaPos))
118 : return false;
119 :
120 0 : if (elem_dim == dim) // full-dimensional element, two different points
121 : {
122 : // extend the ends of the segment to the sides of the element
123 0 : if (ContainsPoint(elem, ls, aaPos))
124 0 : lambda_min = 0.0;
125 0 : if (ContainsPoint(elem, le, aaPos))
126 0 : lambda_max = 1.0;
127 :
128 : // skip the segment if its line cuts the element before or after the segment
129 0 : if (lambda_min < 0.0 || lambda_max > 1.0)
130 0 : return false;
131 :
132 : // compute the corrected ends
133 : VecScaleAdd(le, 1.0, ls, lambda_max, dir);
134 : VecScaleAdd(ls, 1.0, ls, lambda_min, dir);
135 :
136 : // restrict line segment to scv (so that ls and le are the intersection points with the sides of the scv)
137 : MathVector<dim> p;
138 0 : for (size_t i = 0; i < geo.num_scvf(); i++)
139 : {
140 : const typename TFVGeom::SCVF& scvf = geo.scvf(i);
141 : number d_min, d_max, lambda;
142 0 : if (scvf.from() == co)
143 : {
144 : VecSubtract(p, ls, scvf.global_ip());
145 : d_min = VecDot(p, scvf.normal());
146 : VecSubtract(p, le, scvf.global_ip());
147 : d_max = VecDot(p, scvf.normal());
148 0 : if (d_min*d_max < 0.0)
149 : {
150 0 : lambda = fabs(d_min)/(fabs(d_min)+fabs(d_max));
151 0 : VecScaleAdd(p, 1.0-lambda, ls, lambda, le);
152 0 : if (d_max > 0.0)
153 : le = p;
154 : else
155 : ls = p;
156 : }
157 0 : else if (d_min > 0.0) // if so than d_max >= 0, too
158 : return false;
159 : }
160 0 : else if (scvf.to() == co)
161 : {
162 : VecSubtract(p, ls, scvf.global_ip());
163 : d_min = VecDot(p, scvf.normal());
164 : VecSubtract(p, le, scvf.global_ip());
165 : d_max = VecDot(p, scvf.normal());
166 0 : if (d_min*d_max < 0.0)
167 : {
168 0 : lambda = fabs(d_min)/(fabs(d_min)+fabs(d_max));
169 0 : VecScaleAdd(p, 1.0-lambda, ls, lambda, le);
170 0 : if (d_max < 0.0)
171 : le = p;
172 : else
173 : ls = p;
174 : }
175 0 : else if (d_min < 0.0) // if so than d_max <= 0, too
176 : return false;
177 : }
178 : }
179 : }
180 0 : else if (elem_dim + 1 == dim) // low-dimensional element, we accept only one point
181 : {
182 : // skip the segment if its line cuts the element before or after the segment
183 0 : if (lambda_min < 0.0 || lambda_max > 1.0)
184 : return false;
185 :
186 : // exclude the case the the line in the plane
187 0 : if (std::fabs (lambda_max - lambda_min) > SMALL)
188 0 : UG_THROW ("FVLineSourceOrSink: Line source or sink lyes in a low-dimensional element!");
189 :
190 : // check if the point is in the scv
191 : VecScaleAdd(ls, 1.0, point1, lambda_min, dir);
192 0 : for (size_t i = 0; i < geo.num_scvf(); i++)
193 : {
194 : const typename TFVGeom::SCVF& scvf = geo.scvf(i);
195 : MathVector<dim> d;
196 : VecSubtract(d, ls, scvf.global_ip());
197 0 : if (scvf.from() == co)
198 : {
199 0 : if (VecDot(d, scvf.normal()) > 0.0)
200 0 : return false;
201 : }
202 0 : else if (scvf.to() == co)
203 : {
204 0 : if (VecDot(d, scvf.normal()) < 0.0)
205 : return false;
206 : }
207 : }
208 : le = ls;
209 : }
210 : else
211 0 : UG_THROW ("FVLineSourceOrSink: Dimension mismatch with a grid element!");
212 :
213 0 : return true;
214 : };
215 :
216 : /**
217 : * Finds the proper point source or sink starting with a given index
218 : */
219 : template <int dim, typename TPointData, typename TLineData>
220 : template <typename TElem, typename TAAPos, typename TFVGeom>
221 0 : void FVSingularSourcesAndSinks<dim, TPointData, TLineData>::point_iterator<TElem, TAAPos, TFVGeom>::next_sss
222 : (
223 : size_t index ///< index to start the search with
224 : )
225 : {
226 0 : for (m_index = index; m_index < m_sss->num_points (); m_index++)
227 : {
228 : point_sss_type * pnt_sss = m_sss->ListP[m_index].get();
229 0 : if (! m_elem_bbox.contains_point (pnt_sss->position ())) // a quick test
230 0 : continue;
231 0 : if (pnt_sss->corresponds_to (m_elem, m_grid, m_aaPos, m_geo, m_co))
232 : break;
233 : }
234 0 : };
235 :
236 : /**
237 : * Finds the proper point source or sink starting with a given index
238 : */
239 : template <int dim, typename TPointData, typename TLineData>
240 : template <typename TElem, typename TAAPos, typename TFVGeom>
241 0 : void FVSingularSourcesAndSinks<dim, TPointData, TLineData>::line_iterator<TElem, TAAPos, TFVGeom>::next_sss
242 : (
243 : size_t index, ///< index to start the search with
244 : MathVector<dim>& ls, ///< the 1st of the intersection points with the bnd of the scv
245 : MathVector<dim>& le ///< the 2nd of the intersection points with the bnd of the scv
246 : )
247 : {
248 0 : for (m_index = index; m_index < m_sss->num_lines (); m_index++)
249 : {
250 : line_sss_type * line_sss = m_sss->ListL[m_index].get();
251 0 : if (! m_elem_bbox.overlaps_line (line_sss->from_position (), line_sss->to_position ())) // a quick test
252 0 : continue;
253 0 : if (line_sss->corresponds_to (m_elem, m_grid, m_aaPos, m_geo, m_co, ls, le))
254 : break;
255 : }
256 0 : };
257 :
258 : } // end of namespace ug
259 :
260 : /* End of File */
|
