KratosMultiphysics
KRATOS Multiphysics (Kratos) is a framework for building parallel, multi-disciplinary simulation software, aiming at modularity, extensibility, and high performance. Kratos is written in C++, and counts with an extensive Python interface.
edge_data_c2c.h
Go to the documentation of this file.
1 /*
2 ==============================================================================
3 KratosPFEMApplication
4 A library based on:
5 Kratos
6 A General Purpose Software for Multi-Physics Finite Element Analysis
7 Version 1.0 (Released on march 05, 2007).
8 
9 Copyright 2007
10 Pooyan Dadvand, Riccardo Rossi
11 pooyan@cimne.upc.edu
12 rrossi@cimne.upc.edu
13 - CIMNE (International Center for Numerical Methods in Engineering),
14 Gran Capita' s/n, 08034 Barcelona, Spain
15 
16 
17 Permission is hereby granted, free of charge, to any person obtaining
18 a copy of this software and associated documentation files (the
19 "Software"), to deal in the Software without restriction, including
20 without limitation the rights to use, copy, modify, merge, publish,
21 distribute, sublicense and/or sell copies of the Software, and to
22 permit persons to whom the Software is furnished to do so, subject to
23 the following condition:
24 
25 Distribution of this code for any commercial purpose is permissible
26 ONLY BY DIRECT ARRANGEMENT WITH THE COPYRIGHT OWNERS.
27 
28 The above copyright notice and this permission notice shall be
29 included in all copies or substantial portions of the Software.
30 
31 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
32 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
33 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
34 IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
35 CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
36 TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
37 SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
38 
39 ==============================================================================
40  */
41 
42 //
43 // Project Name: Kratos
44 // Last Modified by: $Author: antonia $
45 // Date: $Date: 2009-01-14 08:26:51 $
46 // Revision: $Revision: 1.11 $
47 //
48 //
49 
50 #if !defined(KRATOS_EDGE_DATA_C2C_H_INCLUDED)
51 #define KRATOS_EDGE_DATA_C2C_H_INCLUDED
52 
53 // we suggest defining the following macro
54 #define USE_CONSERVATIVE_FORM_FOR_SCALAR_CONVECTION
55 
56 // we suggest defining the following macro*/*/
57 // #define USE_CONSERVATIVE_FORM_FOR_VECTOR_CONVECTION
58 
59 // System includes
60 #include <string>
61 #include <iostream>
62 #include <algorithm>
63 
64 // External includes
65 
66 // Project includes
67 #include "includes/define.h"
68 #include "includes/model_part.h"
69 #include "includes/node.h"
70 #include "utilities/geometry_utilities.h"
71 #include "free_surface_application.h"
72 #include "utilities/openmp_utils.h"
73 
74 namespace Kratos
75 {
76  // structure definition for fast access to edge data using CSR format
77 
78  template <unsigned int TDim>
79  class EdgesStructureTypeC2C
80  {
81  public:
82  // component ij of the consistent mass matrix (M = Ni * Nj * dOmega)
83  double Mass = 0.0;
84  // components kl of the laplacian matrix of edge ij (L = dNi/dxk * dNj/dxl * dOmega)
85  // double Laplacian;
86  boost::numeric::ublas::bounded_matrix<double, TDim, TDim> LaplacianIJ;
87  // components k of the gradient matrix of edge ij (G = Ni * dNj/dxl * dOmega)
88  array_1d<double, TDim> Ni_DNj = ZeroVector(TDim);
89  // components k of the transposed gradient matrix of edge ij (GT = dNi/dxl * Nj * dOmega)
90  // TRANSPOSED GRADIENT
91  array_1d<double, TDim> DNi_Nj = ZeroVector(TDim);
92 
93  //*************************************************************************************
94  //*************************************************************************************
95  // gradient integrated by parts
96  // RHSi += DNi_Nj pj + Aboundary * pext ==> RHS += Ni_DNj p_j - DNi_Nj p_i
97  // ATTENTION: + Aboundary * pext is NOT included!! it should be included "manually"
98 
99  inline void Add_Gp(array_1d<double, TDim> &destination, const double &p_i, const double &p_j)
100  {
101  for (unsigned int comp = 0; comp < TDim; comp++)
102  destination[comp] -= Ni_DNj[comp] * p_j - DNi_Nj[comp] * p_i;
103  }
104 
105  inline void Sub_Gp(array_1d<double, TDim> &destination, const double &p_i, const double &p_j)
106  {
107  for (unsigned int comp = 0; comp < TDim; comp++)
108  destination[comp] += Ni_DNj[comp] * p_j - DNi_Nj[comp] * p_i;
109  }
110 
111  //*************************************************************************************
112  //*************************************************************************************
113  // gradient
114  // RHSi += Ni_DNj[k]*v[k]
115 
116  inline void Add_D_v(double &destination,
117  const array_1d<double, TDim> &v_i,
118  const array_1d<double, TDim> &v_j)
119  {
120  for (unsigned int comp = 0; comp < TDim; comp++)
121  destination += Ni_DNj[comp] * (v_j[comp] - v_i[comp]);
122  }
123 
124  inline void Sub_D_v(double &destination,
125  const array_1d<double, TDim> &v_i,
126  const array_1d<double, TDim> &v_j)
127  {
128  for (unsigned int comp = 0; comp < TDim; comp++)
129  destination -= Ni_DNj[comp] * (v_j[comp] - v_i[comp]);
130  }
131 
132  //*************************************************************************************
133  //*************************************************************************************
134  // gradient
135  // RHSi += Ni_DNj pj
136 
137  inline void Add_grad_p(array_1d<double, TDim> &destination, const double &p_i, const double &p_j)
138  {
139  for (unsigned int comp = 0; comp < TDim; comp++)
140  destination[comp] += Ni_DNj[comp] * (p_j - p_i);
141  }
142 
143  inline void Sub_grad_p(array_1d<double, TDim> &destination, const double &p_i, const double &p_j)
144  {
145  for (unsigned int comp = 0; comp < TDim; comp++)
146  destination[comp] -= Ni_DNj[comp] * (p_j - p_i);
147  }
148 
149  //*************************************************************************************
150  //*************************************************************************************
151  // gradient
152  // RHSi += DNi_Nj[k]*v[k]
153 
154  inline void Add_div_v(double &destination,
155  const array_1d<double, TDim> &v_i,
156  const array_1d<double, TDim> &v_j)
157  {
158  for (unsigned int comp = 0; comp < TDim; comp++)
159  destination -= Ni_DNj[comp] * v_j[comp] - DNi_Nj[comp] * v_i[comp];
160  }
161 
162  inline void Sub_div_v(double &destination,
163  const array_1d<double, TDim> &v_i,
164  const array_1d<double, TDim> &v_j)
165  {
166  for (unsigned int comp = 0; comp < TDim; comp++)
167  destination += Ni_DNj[comp] * v_j[comp] - DNi_Nj[comp] * v_i[comp];
168  }
169 
170  //*************************************************************************************
171  //*************************************************************************************
172  // gets the trace of the laplacian matrix
173 
174  inline void CalculateScalarLaplacian(double &l_ij)
175  {
176  l_ij = LaplacianIJ(0, 0);
177  for (unsigned int comp = 1; comp < TDim; comp++)
178  l_ij += LaplacianIJ(comp, comp);
179  }
180 
181  inline void Add_ConvectiveContribution(array_1d<double, TDim> &destination,
182  const array_1d<double, TDim> &a_i, const array_1d<double, TDim> &U_i,
183  const array_1d<double, TDim> &a_j, const array_1d<double, TDim> &U_j)
184  {
185  double second = a_i[0] * DNi_Nj[0];
186  double first = a_j[0] * Ni_DNj[0];
187  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
188  {
189  second += a_i[k_comp] * DNi_Nj[k_comp];
190  first += a_j[k_comp] * Ni_DNj[k_comp];
191  }
192  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
193  destination[l_comp] += first * U_j[l_comp] - second * U_i[l_comp];
194  }
195 
196  inline void Sub_ConvectiveContribution(array_1d<double, TDim> &destination,
197  const array_1d<double, TDim> &a_i, const array_1d<double, TDim> &U_i,
198  const array_1d<double, TDim> &a_j, const array_1d<double, TDim> &U_j)
199  {
200  double second = a_i[0] * DNi_Nj[0];
201  double first = a_j[0] * Ni_DNj[0];
202  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
203  {
204  second += a_i[k_comp] * DNi_Nj[k_comp];
205  first += a_j[k_comp] * Ni_DNj[k_comp];
206  }
207  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
208  destination[l_comp] -= first * U_j[l_comp] - second * U_i[l_comp];
209  }
210 
211  inline void Sub_ConvectiveContribution(double &destination,
212  const array_1d<double, TDim> &a_i, const double &phi_i,
213  const array_1d<double, TDim> &a_j, const double &phi_j)
214  {
215 
216 #ifdef USE_CONSERVATIVE_FORM_FOR_SCALAR_CONVECTION
217  double temp = a_i[0] * Ni_DNj[0];
218  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
219  temp += a_i[k_comp] * Ni_DNj[k_comp];
220 
221  destination -= temp * (phi_j - phi_i);
222 #else
223  double aux_i = a_i[0] * Ni_DNj[0];
224  double aux_j = a_j[0] * Ni_DNj[0];
225  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
226  {
227  aux_i += a_i[k_comp] * Ni_DNj[k_comp];
228  aux_j += a_j[k_comp] * Ni_DNj[k_comp];
229  }
230  destination -= aux_j * phi_j - aux_i * phi_i;
231 #endif
232  }
233 
234  inline void Add_ConvectiveContribution(double &destination,
235  const array_1d<double, TDim> &a_i, const double &phi_i,
236  const array_1d<double, TDim> &a_j, const double &phi_j)
237  {
238 
239 #ifdef USE_CONSERVATIVE_FORM_FOR_SCALAR_CONVECTION
240  double temp = a_i[0] * Ni_DNj[0];
241  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
242  temp += a_i[k_comp] * Ni_DNj[k_comp];
243 
244  destination += temp * (phi_j - phi_i);
245 #else
246  double aux_i = a_i[0] * Ni_DNj[0];
247  double aux_j = a_j[0] * Ni_DNj[0];
248  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
249  {
250  aux_i += a_i[k_comp] * Ni_DNj[k_comp];
251  aux_j += a_j[k_comp] * Ni_DNj[k_comp];
252  }
253  destination += aux_j * phi_j - aux_i * phi_i;
254 #endif
255  }
256 
257  //*************************************************************************************
258  //*************************************************************************************
259 
260  inline void CalculateConvectionStabilization_LOW(array_1d<double, TDim> &stab_low,
261  const array_1d<double, TDim> &a_i, const array_1d<double, TDim> &U_i,
262  const array_1d<double, TDim> &a_j, const array_1d<double, TDim> &U_j)
263  {
264  double conv_stab = 0.0;
265  for (unsigned int k_comp = 0; k_comp < TDim; k_comp++)
266  for (unsigned int m_comp = 0; m_comp < TDim; m_comp++)
267  conv_stab += a_i[k_comp] * a_i[m_comp] * LaplacianIJ(k_comp, m_comp);
268  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
269  stab_low[l_comp] = conv_stab * (U_j[l_comp] - U_i[l_comp]);
270 
271  }
272 
273  inline void CalculateConvectionStabilization_LOW(double &stab_low,
274  const array_1d<double, TDim> &a_i, const double &phi_i,
275  const array_1d<double, TDim> &a_j, const double &phi_j)
276  {
277  double conv_stab = 0.0;
278  for (unsigned int k_comp = 0; k_comp < TDim; k_comp++)
279  for (unsigned int m_comp = 0; m_comp < TDim; m_comp++)
280  conv_stab += a_i[k_comp] * a_i[m_comp] * LaplacianIJ(k_comp, m_comp);
281  stab_low = conv_stab * (phi_j - phi_i);
282  }
283  //*************************************************************************************
284  //*************************************************************************************
285 
286  inline void CalculateConvectionStabilization_HIGH(array_1d<double, TDim> &stab_high,
287  const array_1d<double, TDim> &a_i, const array_1d<double, TDim> &pi_i,
288  const array_1d<double, TDim> &a_j, const array_1d<double, TDim> &pi_j)
289  {
290 
291 #ifdef USE_CONSERVATIVE_FORM_FOR_VECTOR_CONVECTION
292  double temp = 0.0;
293  for (unsigned int k_comp = 0; k_comp < TDim; k_comp++)
294  temp += a_i[k_comp] * Ni_DNj[k_comp];
295  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
296  stab_high[l_comp] = -temp * (pi_j[l_comp] - pi_i[l_comp]); // check if the minus sign is correct
297 #else
298  double aux_i = a_i[0] * Ni_DNj[0];
299  double aux_j = a_j[0] * Ni_DNj[0];
300  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
301  {
302  aux_i += a_i[k_comp] * Ni_DNj[k_comp];
303  aux_j += a_j[k_comp] * Ni_DNj[k_comp];
304  }
305  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
306  stab_high[l_comp] = -(aux_j * pi_j[l_comp] - aux_i * pi_i[l_comp]);
307 #endif
308  }
309 
310  inline void CalculateConvectionStabilization_HIGH(double &stab_high,
311  const array_1d<double, TDim> &a_i, const double &pi_i,
312  const array_1d<double, TDim> &a_j, const double &pi_j)
313  {
314 
315 #ifdef USE_CONSERVATIVE_FORM_FOR_SCALAR_CONVECTION
316  double temp = 0.0;
317  for (unsigned int k_comp = 0; k_comp < TDim; k_comp++)
318  temp += a_i[k_comp] * Ni_DNj[k_comp];
319 
320  stab_high = -temp * (pi_j - pi_i); // check if the minus sign is correct
321 #else
322  double aux_i = a_i[0] * Ni_DNj[0];
323  double aux_j = a_j[0] * Ni_DNj[0];
324  for (unsigned int k_comp = 1; k_comp < TDim; k_comp++)
325  {
326  aux_i += a_i[k_comp] * Ni_DNj[k_comp];
327  aux_j += a_j[k_comp] * Ni_DNj[k_comp];
328  }
329 
330  stab_high = -(aux_j * pi_j - aux_i * pi_i);
331 #endif
332  }
333  //*************************************************************************************
334  //*************************************************************************************
335 
336  inline void Add_StabContribution(array_1d<double, TDim> &destination,
337  const double tau, const double beta,
338  const array_1d<double, TDim> &stab_low, const array_1d<double, TDim> &stab_high)
339  {
340  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
341  destination[l_comp] += tau * (stab_low[l_comp] - beta * stab_high[l_comp]);
342  }
343 
344  inline void Add_StabContribution(double &destination,
345  const double tau, const double beta,
346  const double &stab_low, const double &stab_high)
347  {
348  destination += tau * (stab_low - beta * stab_high);
349  }
350 
351  inline void Sub_StabContribution(array_1d<double, TDim> &destination,
352  const double tau, const double beta,
353  const array_1d<double, TDim> &stab_low, const array_1d<double, TDim> &stab_high)
354  {
355  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
356  destination[l_comp] -= tau * (stab_low[l_comp] - beta * stab_high[l_comp]);
357  }
358 
359  inline void Sub_StabContribution(double &destination,
360  const double tau, const double beta,
361  const double &stab_low, const double &stab_high)
362  {
363  destination -= tau * (stab_low - beta * stab_high);
364  }
365 
366  //*************************************************************************************
367  //*************************************************************************************
368 
369  inline void Add_ViscousContribution(array_1d<double, TDim> &destination,
370  const array_1d<double, TDim> &U_i, const double &nu_i,
371  const array_1d<double, TDim> &U_j, const double &nu_j)
372  {
373  // calculate scalar laplacian
374  double L = 0.0;
375  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
376  L += LaplacianIJ(l_comp, l_comp);
377 
378  // double nu_avg = 0.5*(nu_i+nu_j);
379  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
380  destination[l_comp] += nu_i * L * (U_j[l_comp] - U_i[l_comp]);
381  }
382 
383  inline void Sub_ViscousContribution(array_1d<double, TDim> &destination,
384  const array_1d<double, TDim> &U_i, const double &nu_i,
385  const array_1d<double, TDim> &U_j, const double &nu_j)
386  {
387  // calculate scalar laplacian
388  double L = 0.0;
389  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
390  L += LaplacianIJ(l_comp, l_comp);
391 
392  // double nu_avg = 0.5*(nu_i+nu_j);
393  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
394  destination[l_comp] -= nu_i * L * (U_j[l_comp] - U_i[l_comp]);
395  }
396  };
397 
398  // class definition of matrices using CSR format
399 
400  template <unsigned int TDim, class TSparseSpace>
401  class MatrixContainerC2C
402  {
403  public:
404  // name for the self defined structure
405  typedef EdgesStructureTypeC2C<TDim> CSR_Tuple;
406  typedef vector<CSR_Tuple> EdgesVectorType;
407  // name for row start and column index vectors
408  typedef vector<unsigned int> IndicesVectorType;
409  // names for separately stored node based values
410  typedef vector<double> ValuesVectorType;
411  typedef vector<array_1d<double, TDim>> CalcVectorType;
412 
413  // constructor and destructor
414 
415  MatrixContainerC2C(){};
416 
417  ~MatrixContainerC2C(){};
418 
419  // functions to return private values
420 
421  inline unsigned int GetNumberEdges()
422  {
423  return mNumberEdges;
424  }
425 
426  inline EdgesVectorType &GetEdgeValues()
427  {
428  return mNonzeroEdgeValues;
429  }
430 
431  inline IndicesVectorType &GetColumnIndex()
432  {
433  return mColumnIndex;
434  }
435 
436  inline IndicesVectorType &GetRowStartIndex()
437  {
438  return mRowStartIndex;
439  }
440 
441  inline ValuesVectorType &GetLumpedMass()
442  {
443  return mLumpedMassMatrix;
444  }
445 
446  inline ValuesVectorType &GetInvertedMass()
447  {
448  return mInvertedMassMatrix;
449  }
450 
451  inline CalcVectorType &GetDiagGradient()
452  {
453  return mDiagGradientMatrix;
454  }
455 
456  inline ValuesVectorType &GetHmin()
457  {
458  return mHmin;
459  }
460 
461  //********************************************************
462  // function to size and initialize the vector of CSR tuples
463 
464  void ConstructCSRVector(ModelPart &model_part)
465  {
466  KRATOS_TRY
467 
468  // SIZE OF CSR VECTOR
469 
470  // defining the number of nodes and edges
471  int n_nodes = model_part.Nodes().size();
472  // remark: no colouring algorithm is used here (symmetry is neglected)
473  // respectively edge ij is considered different from edge ji
474  mNumberEdges = 0;
475  // counter to assign and get global nodal index
476  int i_node = 0;
477 
478  // counting the edges connecting the nodes
479  for (typename ModelPart::NodesContainerType::iterator node_it = model_part.NodesBegin(); node_it != model_part.NodesEnd(); node_it++)
480  {
481  // counting neighbours of each node
482  mNumberEdges += (node_it->GetValue(NEIGHBOUR_NODES)).size();
483 
484  // assigning global index to each node
485  node_it->FastGetSolutionStepValue(AUX_INDEX) = static_cast<double>(i_node++);
486  }
487  // error message in case number of nodes does not coincide with number of indices
488  if (i_node != n_nodes)
489  KRATOS_WATCH("ERROR - Highest nodal index doesn't coincide with number of nodes!");
490 
491  // allocating memory for block of CSR data - setting to zero for first-touch OpenMP allocation
492  mNonzeroEdgeValues.resize(mNumberEdges);
493  SetToZero(mNonzeroEdgeValues);
494  mColumnIndex.resize(mNumberEdges);
495  SetToZero(mColumnIndex);
496  mRowStartIndex.resize(n_nodes + 1);
497  SetToZero(mRowStartIndex);
498  mLumpedMassMatrix.resize(n_nodes);
499  SetToZero(mLumpedMassMatrix);
500  mInvertedMassMatrix.resize(n_nodes);
501  SetToZero(mInvertedMassMatrix);
502  mDiagGradientMatrix.resize(n_nodes);
503  SetToZero(mDiagGradientMatrix);
504  mHmin.resize(n_nodes);
505  SetToZero(mHmin);
506 
507  // INITIALIZING OF THE CSR VECTOR
508 
509  // temporary variable as the row start index of a node depends on the number of neighbours of the previous one
510  unsigned int row_start_temp = 0;
511 
512  int number_of_threads = ParallelUtilities::GetNumThreads();
513  std::vector<int> row_partition(number_of_threads);
514  OpenMPUtils::DivideInPartitions(model_part.Nodes().size(), number_of_threads, row_partition);
515 
516  for (int k = 0; k < number_of_threads; k++)
517  {
518 #pragma omp parallel
519  if (OpenMPUtils::ThisThread() == k)
520  {
521  // main loop over all nodes
522  for (unsigned int aux_i = static_cast<unsigned int>(row_partition[k]); aux_i < static_cast<unsigned int>(row_partition[k + 1]); aux_i++)
523  {
524  typename ModelPart::NodesContainerType::iterator node_it = model_part.NodesBegin() + aux_i;
525 
526  // getting the global index of the node
527  i_node = static_cast<unsigned int>(node_it->FastGetSolutionStepValue(AUX_INDEX));
528 
529  // determining its neighbours
530  GlobalPointersVector<Node> &neighb_nodes = node_it->GetValue(NEIGHBOUR_NODES);
531 
532  // number of neighbours of node i determines row start index for the following node
533  unsigned int n_neighbours = neighb_nodes.size();
534  // DIAGONAL TERMS
535 
536  // reserving memory for work array
537  std::vector<unsigned int> work_array;
538  work_array.reserve(n_neighbours);
539  // DIAGONAL TERMS
540 
541  // nested loop over the neighbouring nodes
542  for (GlobalPointersVector<Node>::iterator neighb_it = neighb_nodes.begin(); neighb_it != neighb_nodes.end(); neighb_it++)
543  {
544  // getting global index of the neighbouring node
545  work_array.push_back(static_cast<unsigned int>(neighb_it->FastGetSolutionStepValue(AUX_INDEX)));
546  }
547  // reordering neighbours following their global indices
548  std::sort(work_array.begin(), work_array.end());
549 
550  // setting current row start index
551  mRowStartIndex[i_node] = row_start_temp;
552  // nested loop over the by now ordered neighbours
553  for (unsigned int counter = 0; counter < n_neighbours; counter++)
554  {
555  // getting global index of the neighbouring node
556  unsigned int j_neighbour = work_array[counter];
557 
558  // calculating CSR index
559  unsigned int csr_index = mRowStartIndex[i_node] + counter;
560 
561  // saving column index j of the original matrix
562  mColumnIndex[csr_index] = j_neighbour;
563  }
564  // preparing row start index for next node
565  row_start_temp += n_neighbours;
566  }
567  }
568  }
569  // adding last entry (necessary for abort criterion of loops)
570  mRowStartIndex[n_nodes] = mNumberEdges;
571 
572  // INITIALIZING NODE BASED VALUES
573  // set the heights to a huge number
574  IndexPartition<unsigned int>(n_nodes).for_each([&](unsigned int i_node){
575  mHmin[i_node] = 1e10;
576  });
577 
578  KRATOS_CATCH("")
579  }
580 
581  //*********************************
582  // function to precalculate CSR data
583 
584  void BuildCSRData(ModelPart &model_part)
585  {
586  KRATOS_TRY
587 
588  // PRECALCULATING CSR DATA
589 
590  // defining temporary local variables for elementwise addition
591  // shape functions
592  array_1d<double, TDim + 1> N;
593  // shape function derivatives
594  boost::numeric::ublas::bounded_matrix<double, TDim + 1, TDim> dN_dx;
595  // volume
596  double volume;
597  // weighting factor
598  double weighting_factor = 1.0 / static_cast<double>(TDim + 1);
599  // elemental matrices
600  boost::numeric::ublas::bounded_matrix<double, TDim + 1, TDim + 1> mass_consistent;
601  // boost::numeric::ublas::bounded_matrix <double, TDim+1,TDim+1> laplacian;
602  array_1d<double, TDim + 1> mass_lumped;
603  // global indices of elemental nodes
604  array_1d<unsigned int, TDim + 1> nodal_indices;
605 
606  array_1d<double, TDim + 1> heights;
607 
608  // loop over all elements
609  for (typename ModelPart::ElementsContainerType::iterator elem_it = model_part.ElementsBegin(); elem_it != model_part.ElementsEnd(); elem_it++)
610  {
611  // LOCAL ELEMENTWISE CALCULATIONS
612 
613  // getting geometry data of the element
614  GeometryUtils::CalculateGeometryData(elem_it->GetGeometry(), dN_dx, N, volume);
615 
616  // calculate lenght of the heights of the element
617  for (unsigned int ie_node = 0; ie_node <= TDim; ie_node++)
618  {
619  heights[ie_node] = dN_dx(ie_node, 0) * dN_dx(ie_node, 0);
620  for (unsigned int comp = 1; comp < TDim; comp++)
621  {
622  heights[ie_node] += dN_dx(ie_node, comp) * dN_dx(ie_node, comp);
623  }
624  heights[ie_node] = 1.0 / sqrt(heights[ie_node]);
625  }
626 
627  // setting up elemental mass matrices
628  CalculateMassMatrix(mass_consistent, volume);
629  noalias(mass_lumped) = ZeroVector(TDim + 1);
630  for (unsigned int ie_node = 0; ie_node <= TDim; ie_node++)
631  {
632  for (unsigned int je_node = 0; je_node <= TDim; je_node++)
633  {
634  // mass_consistent(ie_node,je_node) = N(ie_node) * N(je_node) * volume;
635  mass_lumped[ie_node] += mass_consistent(ie_node, je_node);
636  }
637  // mass_lumped[ie_node] = volume * N[ie_node];
638  }
639 
640  double weighted_volume = volume * weighting_factor;
641 
642  // ASSEMBLING GLOBAL DATA STRUCTURE
643 
644  // loop over the nodes of the element to determine their global indices
645  for (unsigned int ie_node = 0; ie_node <= TDim; ie_node++)
646  nodal_indices[ie_node] = static_cast<unsigned int>(elem_it->GetGeometry()[ie_node].FastGetSolutionStepValue(AUX_INDEX));
647 
648  // assembling global "edge matrices" by adding local contributions
649  for (unsigned int ie_node = 0; ie_node <= TDim; ie_node++)
650  {
651  // check the heights and change the value if minimal is found
652  if (mHmin[nodal_indices[ie_node]] > heights[ie_node])
653  mHmin[nodal_indices[ie_node]] = heights[ie_node];
654 
655  for (unsigned int je_node = 0; je_node <= TDim; je_node++)
656  {
657  // remark: there is no edge linking node i with itself!
658  // DIAGONAL TERMS
659  if (ie_node != je_node)
660  {
661  // calculating CSR index from global index
662  unsigned int csr_index = GetCSRIndex(nodal_indices[ie_node], nodal_indices[je_node]);
663 
664  // assigning precalculated element data to the referring edges
665  // contribution to edge mass
666  mNonzeroEdgeValues[csr_index].Mass += mass_consistent(ie_node, je_node);
667 
668  // contribution to edge laplacian
669  boost::numeric::ublas::bounded_matrix<double, TDim, TDim> &laplacian = mNonzeroEdgeValues[csr_index].LaplacianIJ;
670  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
671  for (unsigned int k_comp = 0; k_comp < TDim; k_comp++)
672  laplacian(l_comp, k_comp) += dN_dx(ie_node, l_comp) * dN_dx(je_node, k_comp) * volume;
673 
674  // contribution to edge gradient
675  array_1d<double, TDim> &gradient = mNonzeroEdgeValues[csr_index].Ni_DNj;
676  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
677  // gradient[l_comp] += dN_dx(je_node,l_comp);
678  gradient[l_comp] += dN_dx(je_node, l_comp) * weighted_volume;
679  // TRANSPOSED GRADIENT
680  // contribution to transposed edge gradient
681  array_1d<double, TDim> &transp_gradient = mNonzeroEdgeValues[csr_index].DNi_Nj;
682  for (unsigned int l_comp = 0; l_comp < TDim; l_comp++)
683  transp_gradient[l_comp] += dN_dx(ie_node, l_comp) * weighted_volume;
684  }
685  }
686  }
687 
688  // assembling node based vectors
689  for (unsigned int ie_node = 0; ie_node <= TDim; ie_node++)
690  // diagonal of the global lumped mass matrix
691  mLumpedMassMatrix[nodal_indices[ie_node]] += mass_lumped[ie_node];
692  for (unsigned int ie_node = 0; ie_node <= TDim; ie_node++)
693  {
694  // diagonal of the global gradient matrix
695  array_1d<double, TDim> &gradient = mDiagGradientMatrix[nodal_indices[ie_node]];
696  for (unsigned int component = 0; component < TDim; component++)
697  gradient[component] += dN_dx(ie_node, component) * weighted_volume;
698  }
699  }
700 
701  // copy mass matrix to inverted mass matrix
702  for (unsigned int inode = 0; inode < mLumpedMassMatrix.size(); inode++)
703  {
704  mInvertedMassMatrix[inode] = mLumpedMassMatrix[inode];
705  }
706 
707  // perform MPI syncronization between the domains
708 
709  // calculating inverted mass matrix (this requires syncronization for MPI paraellelism
710  for (unsigned int inode = 0; inode < mInvertedMassMatrix.size(); inode++)
711  {
712  mInvertedMassMatrix[inode] = 1.0 / mInvertedMassMatrix[inode];
713  }
714 
715  KRATOS_CATCH("")
716  }
717 
718  //******************************************
719  // function to calculate CSR index of edge ij
720 
721  unsigned int GetCSRIndex(unsigned int NodeI, unsigned int NeighbourJ)
722  {
723  KRATOS_TRY
724 
725  // index indicating data position of edge ij
726  unsigned int csr_index;
727  // searching for coincidence of stored column index and neighbour index j
728  for (csr_index = mRowStartIndex[NodeI]; csr_index != mRowStartIndex[NodeI + 1]; csr_index++)
729  if (mColumnIndex[csr_index] == NeighbourJ)
730  break;
731 
732  // returning CSR index of edge ij
733  return csr_index;
734 
735  KRATOS_CATCH("")
736  }
737 
738  //***********************************************
739  // function to get pointer to CSR tuple of edge ij
740 
741  CSR_Tuple *GetTuplePointer(unsigned int NodeI, unsigned int NeighbourJ)
742  {
743  KRATOS_TRY
744 
745  // index indicating data position of edge ij
746  unsigned int csr_index;
747  // searching for coincidence of stored column index and neighbour index j
748  for (csr_index = mRowStartIndex[NodeI]; csr_index != mRowStartIndex[NodeI + 1]; csr_index++)
749  if (mColumnIndex[csr_index] == NeighbourJ)
750  break;
751 
752  // returning pointer to CSR tuple of edge ij
753  return &mNonzeroEdgeValues[csr_index];
754 
755  KRATOS_CATCH("")
756  }
757 
758  //*******************************
759  // function to free dynamic memory
760 
761  void Clear()
762  {
763  KRATOS_TRY
764 
765  mNonzeroEdgeValues.clear();
766  mColumnIndex.clear();
767  mRowStartIndex.clear();
768  mInvertedMassMatrix.clear();
769  mLumpedMassMatrix.clear();
770  mDiagGradientMatrix.clear();
771  mHmin.clear();
772 
773  KRATOS_CATCH("")
774  }
775  //****************************
776  // functions to access database
777  //(note that this is already thought for parallel;
778  // for a single processor this could be done in a faster way)
779 
780  void FillCoordinatesFromDatabase(CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
781  {
782 
783  KRATOS_TRY
784 
785  // loop over all nodes
786  int n_nodes = rNodes.size();
787  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
788 
789 #pragma omp parallel for firstprivate(n_nodes, it_begin)
790  for (int i = 0; i < n_nodes; i++)
791  {
792  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
793 
794  // get the global index of node i
795  unsigned int i_node = i;
796 
797  // save value in the destination vector
798  for (unsigned int component = 0; component < TDim; component++)
799  (rDestination[i_node])[component] = (*node_it)[component];
800  }
801 
802  KRATOS_CATCH("");
803  }
804 
805  //****************************
806  // functions to access database
807  //(note that this is already thought for parallel;
808  // for a single processor this could be done in a faster way)
809 
810  void FillVectorFromDatabase(Variable<array_1d<double, 3>> &rVariable, CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
811  {
812 
813  KRATOS_TRY
814 
815  // loop over all nodes
816 
817  int n_nodes = rNodes.size();
818 
819  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
820 
821  unsigned int var_pos = it_begin->pGetVariablesList()->Index(rVariable);
822 
823 #pragma omp parallel for firstprivate(n_nodes, it_begin, var_pos)
824  for (int i = 0; i < n_nodes; i++)
825  {
826  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
827 
828  // get the global index of node i
829  unsigned int i_node = i;
830 
831  // get the requested value in vector form
832  array_1d<double, 3> &vector = node_it->FastGetCurrentSolutionStepValue(rVariable, var_pos);
833  // save value in the destination vector
834  for (unsigned int component = 0; component < TDim; component++)
835  (rDestination[i_node])[component] = vector[component];
836  }
837 
838  KRATOS_CATCH("");
839  }
840 
841  void FillOldVectorFromDatabase(Variable<array_1d<double, 3>> &rVariable, CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
842  {
843 
844  KRATOS_TRY
845 
846  // loop over all nodes
847  int n_nodes = rNodes.size();
848 
849  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
850 
851  unsigned int var_pos = it_begin->pGetVariablesList()->Index(rVariable);
852 
853 #pragma omp parallel for firstprivate(n_nodes, it_begin, var_pos)
854  for (int i = 0; i < n_nodes; i++)
855  {
856  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
857 
858  // get the global index of node i
859  unsigned int i_node = i;
860 
861  // get the requested value in vector form
862  array_1d<double, 3> &vector = node_it->FastGetSolutionStepValue(rVariable, 1, var_pos);
863  // save value in the destination vector
864  for (unsigned int component = 0; component < TDim; component++)
865  (rDestination[i_node])[component] = vector[component];
866  }
867 
868  KRATOS_CATCH("");
869  }
870 
871  void FillScalarFromDatabase(Variable<double> &rVariable, ValuesVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
872  {
873  KRATOS_TRY
874 
875  // loop over all nodes
876  int n_nodes = rNodes.size();
877 
878  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
879 
880  unsigned int var_pos = it_begin->pGetVariablesList()->Index(rVariable);
881 
882 #pragma omp parallel for firstprivate(n_nodes, it_begin, var_pos)
883  for (int i = 0; i < n_nodes; i++)
884  {
885  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
886 
887  // get the global index of node i
888  unsigned int i_node = i;
889 
890  // get the requested scalar value
891  double &scalar = node_it->FastGetCurrentSolutionStepValue(rVariable, var_pos);
892  // save value in the destination vector
893  rDestination[i_node] = scalar;
894  }
895 
896  KRATOS_CATCH("");
897  }
898 
899  void FillOldScalarFromDatabase(Variable<double> &rVariable, ValuesVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
900  {
901  KRATOS_TRY
902 
903  int n_nodes = rNodes.size();
904  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
905 
906  unsigned int var_pos = it_begin->pGetVariablesList()->Index(rVariable);
907 
908 #pragma omp parallel for firstprivate(n_nodes, it_begin, var_pos)
909  for (int i = 0; i < n_nodes; i++)
910  {
911  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
912 
913  // get the global index of node i
914  // // unsigned int i_node = static_cast<unsigned int>(node_it->FastGetSolutionStepValue(AUX_INDEX));
915  unsigned int i_node = i;
916 
917  // get the requested scalar value
918  double &scalar = node_it->FastGetSolutionStepValue(rVariable, 1, var_pos);
919  // save value in the destination vector
920  rDestination[i_node] = scalar;
921  }
922 
923  KRATOS_CATCH("");
924  }
925 
926  void WriteVectorToDatabase(Variable<array_1d<double, 3>> &rVariable, CalcVectorType &rOrigin, ModelPart::NodesContainerType &rNodes)
927  {
928  KRATOS_TRY
929 
930  // loop over all nodes
931  int n_nodes = rNodes.size();
932  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
933 
934  unsigned int var_pos = it_begin->pGetVariablesList()->Index(rVariable);
935 
936 #pragma omp parallel for firstprivate(n_nodes, it_begin, var_pos)
937  for (int i = 0; i < n_nodes; i++)
938  {
939  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
940 
941  // get the global index of node i
942  unsigned int i_node = i;
943 
944  // get reference of destination
945  array_1d<double, 3> &vector = node_it->FastGetCurrentSolutionStepValue(rVariable, var_pos);
946  // save vector in database
947  for (unsigned int component = 0; component < TDim; component++)
948  vector[component] = (rOrigin[i_node])[component];
949  }
950 
951  KRATOS_CATCH("");
952  }
953 
954  void WriteScalarToDatabase(Variable<double> &rVariable, ValuesVectorType &rOrigin, ModelPart::NodesContainerType &rNodes)
955  {
956  KRATOS_TRY
957 
958  // loop over all nodes
959  int n_nodes = rNodes.size();
960  ModelPart::NodesContainerType::iterator it_begin = rNodes.begin();
961 
962  unsigned int var_pos = it_begin->pGetVariablesList()->Index(rVariable);
963 
964 #pragma omp parallel for firstprivate(n_nodes, it_begin, var_pos)
965  for (int i = 0; i < n_nodes; i++)
966  {
967  ModelPart::NodesContainerType::iterator node_it = it_begin + i;
968 
969  // get the global index of node i
970  int i_node = i;
971 
972  // get reference of destination
973  double &scalar = node_it->FastGetCurrentSolutionStepValue(rVariable, var_pos);
974  // save scalar in database
975  scalar = rOrigin[i_node];
976  }
977 
978  KRATOS_CATCH("");
979  }
980 
981  //*********************************************************************
982 
983  void Add_Minv_value(
984  CalcVectorType &destination,
985  const CalcVectorType &origin1,
986  const double value,
987  const ValuesVectorType &Minv_vec,
988  const CalcVectorType &origin)
989  {
990  KRATOS_TRY
991 
992  int loop_size = destination.size();
993 
994  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
995  array_1d<double, TDim> &dest = destination[i_node];
996  const double m_inv = Minv_vec[i_node];
997  const array_1d<double, TDim> &origin_vec1 = origin1[i_node];
998  const array_1d<double, TDim> &origin_value = origin[i_node];
999 
1000  double temp = value * m_inv;
1001  for (unsigned int comp = 0; comp < TDim; comp++)
1002  dest[comp] = origin_vec1[comp] + temp * origin_value[comp];
1003  });
1004 
1005  KRATOS_CATCH("")
1006  }
1007 
1008  void Add_Minv_value(
1009  ValuesVectorType &destination,
1010  const ValuesVectorType &origin1,
1011  const double value,
1012  const ValuesVectorType &Minv_vec,
1013  const ValuesVectorType &origin)
1014  {
1015  KRATOS_TRY
1016 
1017  int loop_size = destination.size();
1018 
1019  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1020  double &dest = destination[i_node];
1021  const double m_inv = Minv_vec[i_node];
1022  const double &origin_vec1 = origin1[i_node];
1023  const double &origin_value = origin[i_node];
1024 
1025  double temp = value * m_inv;
1026  dest = origin_vec1 + temp * origin_value;
1027  });
1028 
1029  KRATOS_CATCH("")
1030  }
1031 
1032  //**********************************************************************
1033 
1034  void AllocateAndSetToZero(CalcVectorType &data_vector, int size)
1035  {
1036  data_vector.resize(size);
1037  int loop_size = size;
1038 
1039  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1040  array_1d<double, TDim> &aaa = data_vector[i_node];
1041  for (unsigned int comp = 0; comp < TDim; comp++)
1042  aaa[comp] = 0.0;
1043  });
1044  }
1045 
1046  void AllocateAndSetToZero(ValuesVectorType &data_vector, int size)
1047  {
1048  data_vector.resize(size);
1049  int loop_size = size;
1050 
1051  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1052  data_vector[i_node] = 0.0;
1053  });
1054  }
1055 
1056  //**********************************************************************
1057 
1058  void SetToZero(EdgesVectorType &data_vector)
1059  {
1060  int loop_size = data_vector.size();
1061 
1062  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1063  // initializing the CSR vector entries with zero
1064  data_vector[i_node].Mass = 0.0;
1065  noalias(data_vector[i_node].LaplacianIJ) = ZeroMatrix(TDim, TDim);
1066  noalias(data_vector[i_node].Ni_DNj) = ZeroVector(TDim);
1067  noalias(data_vector[i_node].DNi_Nj) = ZeroVector(TDim);
1068  });
1069  }
1070 
1071  void SetToZero(IndicesVectorType &data_vector)
1072  {
1073  int loop_size = data_vector.size();
1074 
1075  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1076  data_vector[i_node] = 0.0;
1077  });
1078  }
1079 
1080  void SetToZero(CalcVectorType &data_vector)
1081  {
1082  int loop_size = data_vector.size();
1083 
1084  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1085  array_1d<double, TDim> &aaa = data_vector[i_node];
1086  for (unsigned int comp = 0; comp < TDim; comp++)
1087  aaa[comp] = 0.0;
1088  });
1089  }
1090 
1091  void SetToZero(ValuesVectorType &data_vector)
1092  {
1093  int loop_size = data_vector.size();
1094 
1095  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1096  data_vector[i_node] = 0.0;
1097  });
1098  }
1099 
1100  //**********************************************************************
1101 
1102  void AssignVectorToVector(const CalcVectorType &origin,
1103  CalcVectorType &destination)
1104  {
1105  int loop_size = origin.size();
1106 
1107  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1108  const array_1d<double, TDim> &orig = origin[i_node];
1109  array_1d<double, TDim> &dest = destination[i_node];
1110  for (unsigned int comp = 0; comp < TDim; comp++)
1111  dest[comp] = orig[comp];
1112  });
1113  }
1114 
1115  void AssignVectorToVector(const ValuesVectorType &origin,
1116  ValuesVectorType &destination)
1117  {
1118  int loop_size = origin.size();
1119 
1120  IndexPartition<unsigned int>(loop_size).for_each([&](unsigned int i_node){
1121  destination[i_node] = origin[i_node];
1122  });
1123  }
1124 
1125  private:
1126  // number of edges
1127  unsigned int mNumberEdges;
1128 
1129  // CSR data vector for storage of the G, L and consistent M components of edge ij
1130  EdgesVectorType mNonzeroEdgeValues;
1131 
1132  // vector to store column indices of nonzero matrix elements for each row
1133  IndicesVectorType mColumnIndex;
1134 
1135  // index vector to access the start of matrix row i in the column vector
1136  IndicesVectorType mRowStartIndex;
1137 
1138  // inverse of the mass matrix ... for parallel calculation each subdomain should contain this correctly calculated (including contributions of the neighbours)
1139  ValuesVectorType mInvertedMassMatrix;
1140 
1141  // minimum height around one node
1142  ValuesVectorType mHmin;
1143 
1144  // lumped mass matrix (separately stored due to lack of diagonal elements of the consistent mass matrix)
1145  ValuesVectorType mLumpedMassMatrix;
1146  // diagonal of the gradient matrix (separately stored due to special calculations)
1147  CalcVectorType mDiagGradientMatrix;
1148 
1149  //*******************************************
1150  // functions to set up elemental mass matrices
1151 
1152  void CalculateMassMatrix(boost::numeric::ublas::bounded_matrix<double, 3, 3> &mass_consistent, double volume)
1153  {
1154  for (unsigned int i_node = 0; i_node <= TDim; i_node++)
1155  {
1156  // diagonal terms
1157  mass_consistent(i_node, i_node) = 0.16666666666666666667 * volume; // 1/6
1158  // non-diagonal terms
1159  double temp = 0.08333333333333333333 * volume; // 1/12
1160  for (unsigned int j_neighbour = i_node + 1; j_neighbour <= TDim; j_neighbour++)
1161  {
1162  // taking advantage of symmetry
1163  mass_consistent(i_node, j_neighbour) = temp;
1164  mass_consistent(j_neighbour, i_node) = temp;
1165  }
1166  }
1167  }
1168 
1169  void CalculateMassMatrix(boost::numeric::ublas::bounded_matrix<double, 4, 4> &mass_consistent, double volume)
1170  {
1171  for (unsigned int i_node = 0; i_node <= TDim; i_node++)
1172  {
1173  // diagonal terms
1174  mass_consistent(i_node, i_node) = 0.1 * volume;
1175  // non-diagonal terms
1176  double temp = 0.05 * volume;
1177  for (unsigned int j_neighbour = i_node + 1; j_neighbour <= TDim; j_neighbour++)
1178  {
1179  // taking advantage of symmetry
1180  mass_consistent(i_node, j_neighbour) = temp;
1181  mass_consistent(j_neighbour, i_node) = temp;
1182  }
1183  }
1184  }
1185  };
1186 
1187 } // namespace Kratos
1188 
1189 #endif // KRATOS_EDGE_DATA_C2C_H_INCLUDED defined
Kratos::EdgesStructureTypeC2C
Definition: edge_data_c2c.h:80
Kratos::EdgesStructureTypeC2C::CalculateScalarLaplacian
void CalculateScalarLaplacian(double &l_ij)
Definition: edge_data_c2c.h:174
Kratos::EdgesStructureTypeC2C::CalculateConvectionStabilization_LOW
void CalculateConvectionStabilization_LOW(array_1d< double, TDim > &stab_low, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &U_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &U_j)
Definition: edge_data_c2c.h:260
Kratos::EdgesStructureTypeC2C::Sub_StabContribution
void Sub_StabContribution(array_1d< double, TDim > &destination, const double tau, const double beta, const array_1d< double, TDim > &stab_low, const array_1d< double, TDim > &stab_high)
Definition: edge_data_c2c.h:351
Kratos::EdgesStructureTypeC2C::Ni_DNj
array_1d< double, TDim > Ni_DNj
Definition: edge_data_c2c.h:88
Kratos::EdgesStructureTypeC2C::Add_div_v
void Add_div_v(double &destination, const array_1d< double, TDim > &v_i, const array_1d< double, TDim > &v_j)
Definition: edge_data_c2c.h:154
Kratos::EdgesStructureTypeC2C::Sub_Gp
void Sub_Gp(array_1d< double, TDim > &destination, const double &p_i, const double &p_j)
Definition: edge_data_c2c.h:105
Kratos::EdgesStructureTypeC2C::Sub_div_v
void Sub_div_v(double &destination, const array_1d< double, TDim > &v_i, const array_1d< double, TDim > &v_j)
Definition: edge_data_c2c.h:162
Kratos::EdgesStructureTypeC2C::CalculateConvectionStabilization_LOW
void CalculateConvectionStabilization_LOW(double &stab_low, const array_1d< double, TDim > &a_i, const double &phi_i, const array_1d< double, TDim > &a_j, const double &phi_j)
Definition: edge_data_c2c.h:273
Kratos::EdgesStructureTypeC2C::Add_StabContribution
void Add_StabContribution(double &destination, const double tau, const double beta, const double &stab_low, const double &stab_high)
Definition: edge_data_c2c.h:344
Kratos::EdgesStructureTypeC2C::CalculateConvectionStabilization_HIGH
void CalculateConvectionStabilization_HIGH(double &stab_high, const array_1d< double, TDim > &a_i, const double &pi_i, const array_1d< double, TDim > &a_j, const double &pi_j)
Definition: edge_data_c2c.h:310
Kratos::EdgesStructureTypeC2C::Add_ConvectiveContribution
void Add_ConvectiveContribution(array_1d< double, TDim > &destination, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &U_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &U_j)
Definition: edge_data_c2c.h:181
Kratos::EdgesStructureTypeC2C::Sub_StabContribution
void Sub_StabContribution(double &destination, const double tau, const double beta, const double &stab_low, const double &stab_high)
Definition: edge_data_c2c.h:359
Kratos::EdgesStructureTypeC2C::Sub_ConvectiveContribution
void Sub_ConvectiveContribution(double &destination, const array_1d< double, TDim > &a_i, const double &phi_i, const array_1d< double, TDim > &a_j, const double &phi_j)
Definition: edge_data_c2c.h:211
Kratos::EdgesStructureTypeC2C::Sub_ConvectiveContribution
void Sub_ConvectiveContribution(array_1d< double, TDim > &destination, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &U_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &U_j)
Definition: edge_data_c2c.h:196
Kratos::EdgesStructureTypeC2C::DNi_Nj
array_1d< double, TDim > DNi_Nj
Definition: edge_data_c2c.h:91
Kratos::EdgesStructureTypeC2C::Add_Gp
void Add_Gp(array_1d< double, TDim > &destination, const double &p_i, const double &p_j)
Definition: edge_data_c2c.h:99
Kratos::EdgesStructureTypeC2C::Add_StabContribution
void Add_StabContribution(array_1d< double, TDim > &destination, const double tau, const double beta, const array_1d< double, TDim > &stab_low, const array_1d< double, TDim > &stab_high)
Definition: edge_data_c2c.h:336
Kratos::EdgesStructureTypeC2C::Sub_grad_p
void Sub_grad_p(array_1d< double, TDim > &destination, const double &p_i, const double &p_j)
Definition: edge_data_c2c.h:143
Kratos::EdgesStructureTypeC2C::Sub_ViscousContribution
void Sub_ViscousContribution(array_1d< double, TDim > &destination, const array_1d< double, TDim > &U_i, const double &nu_i, const array_1d< double, TDim > &U_j, const double &nu_j)
Definition: edge_data_c2c.h:383
Kratos::EdgesStructureTypeC2C::Add_grad_p
void Add_grad_p(array_1d< double, TDim > &destination, const double &p_i, const double &p_j)
Definition: edge_data_c2c.h:137
Kratos::EdgesStructureTypeC2C::Mass
double Mass
Definition: edge_data_c2c.h:83
Kratos::EdgesStructureTypeC2C::CalculateConvectionStabilization_HIGH
void CalculateConvectionStabilization_HIGH(array_1d< double, TDim > &stab_high, const array_1d< double, TDim > &a_i, const array_1d< double, TDim > &pi_i, const array_1d< double, TDim > &a_j, const array_1d< double, TDim > &pi_j)
Definition: edge_data_c2c.h:286
Kratos::EdgesStructureTypeC2C::Sub_D_v
void Sub_D_v(double &destination, const array_1d< double, TDim > &v_i, const array_1d< double, TDim > &v_j)
Definition: edge_data_c2c.h:124
Kratos::EdgesStructureTypeC2C::Add_ViscousContribution
void Add_ViscousContribution(array_1d< double, TDim > &destination, const array_1d< double, TDim > &U_i, const double &nu_i, const array_1d< double, TDim > &U_j, const double &nu_j)
Definition: edge_data_c2c.h:369
Kratos::EdgesStructureTypeC2C::Add_ConvectiveContribution
void Add_ConvectiveContribution(double &destination, const array_1d< double, TDim > &a_i, const double &phi_i, const array_1d< double, TDim > &a_j, const double &phi_j)
Definition: edge_data_c2c.h:234
Kratos::EdgesStructureTypeC2C::Add_D_v
void Add_D_v(double &destination, const array_1d< double, TDim > &v_i, const array_1d< double, TDim > &v_j)
Definition: edge_data_c2c.h:116
Kratos::EdgesStructureTypeC2C::LaplacianIJ
boost::numeric::ublas::bounded_matrix< double, TDim, TDim > LaplacianIJ
Definition: edge_data_c2c.h:86
Kratos::GeometryUtils::CalculateGeometryData
static void CalculateGeometryData(const GeometryType &rGeometry, BoundedMatrix< double, 4, 3 > &rDN_DX, array_1d< double, 4 > &rN, double &rVolume)
This function is designed to compute the shape function derivatives, shape functions and volume in 3D...
Definition: geometry_utilities.h:176
Kratos::GlobalPointersVector
This class is a vector which stores global pointers.
Definition: global_pointers_vector.h:61
Kratos::GlobalPointersVector::iterator
boost::indirect_iterator< typename TContainerType::iterator > iterator
Definition: global_pointers_vector.h:79
Kratos::GlobalPointersVector::begin
iterator begin()
Definition: global_pointers_vector.h:221
Kratos::GlobalPointersVector::size
size_type size() const
Definition: global_pointers_vector.h:307
Kratos::GlobalPointersVector::end
iterator end()
Definition: global_pointers_vector.h:229
Kratos::IndexPartition
This class is useful for index iteration over containers.
Definition: parallel_utilities.h:451
Kratos::IndexPartition::for_each
void for_each(TUnaryFunction &&f)
Definition: parallel_utilities.h:514
Kratos::MatrixContainerC2C
Definition: edge_data_c2c.h:402
Kratos::MatrixContainerC2C::GetLumpedMass
ValuesVectorType & GetLumpedMass()
Definition: edge_data_c2c.h:441
Kratos::MatrixContainerC2C::AssignVectorToVector
void AssignVectorToVector(const ValuesVectorType &origin, ValuesVectorType &destination)
Definition: edge_data_c2c.h:1115
Kratos::MatrixContainerC2C::CalcVectorType
vector< array_1d< double, TDim > > CalcVectorType
Definition: edge_data_c2c.h:411
Kratos::MatrixContainerC2C::Add_Minv_value
void Add_Minv_value(CalcVectorType &destination, const CalcVectorType &origin1, const double value, const ValuesVectorType &Minv_vec, const CalcVectorType &origin)
Definition: edge_data_c2c.h:983
Kratos::MatrixContainerC2C::FillCoordinatesFromDatabase
void FillCoordinatesFromDatabase(CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:780
Kratos::MatrixContainerC2C::~MatrixContainerC2C
~MatrixContainerC2C()
Definition: edge_data_c2c.h:417
Kratos::MatrixContainerC2C::FillVectorFromDatabase
void FillVectorFromDatabase(Variable< array_1d< double, 3 >> &rVariable, CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:810
Kratos::MatrixContainerC2C::GetDiagGradient
CalcVectorType & GetDiagGradient()
Definition: edge_data_c2c.h:451
Kratos::MatrixContainerC2C::CSR_Tuple
EdgesStructureTypeC2C< TDim > CSR_Tuple
Definition: edge_data_c2c.h:405
Kratos::MatrixContainerC2C::GetNumberEdges
unsigned int GetNumberEdges()
Definition: edge_data_c2c.h:421
Kratos::MatrixContainerC2C::ValuesVectorType
vector< double > ValuesVectorType
Definition: edge_data_c2c.h:410
Kratos::MatrixContainerC2C::FillOldVectorFromDatabase
void FillOldVectorFromDatabase(Variable< array_1d< double, 3 >> &rVariable, CalcVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:841
Kratos::MatrixContainerC2C::EdgesVectorType
vector< CSR_Tuple > EdgesVectorType
Definition: edge_data_c2c.h:406
Kratos::MatrixContainerC2C::GetEdgeValues
EdgesVectorType & GetEdgeValues()
Definition: edge_data_c2c.h:426
Kratos::MatrixContainerC2C::GetColumnIndex
IndicesVectorType & GetColumnIndex()
Definition: edge_data_c2c.h:431
Kratos::MatrixContainerC2C::FillOldScalarFromDatabase
void FillOldScalarFromDatabase(Variable< double > &rVariable, ValuesVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:899
Kratos::MatrixContainerC2C::GetRowStartIndex
IndicesVectorType & GetRowStartIndex()
Definition: edge_data_c2c.h:436
Kratos::MatrixContainerC2C::Add_Minv_value
void Add_Minv_value(ValuesVectorType &destination, const ValuesVectorType &origin1, const double value, const ValuesVectorType &Minv_vec, const ValuesVectorType &origin)
Definition: edge_data_c2c.h:1008
Kratos::MatrixContainerC2C::GetHmin
ValuesVectorType & GetHmin()
Definition: edge_data_c2c.h:456
Kratos::MatrixContainerC2C::SetToZero
void SetToZero(EdgesVectorType &data_vector)
Definition: edge_data_c2c.h:1058
Kratos::MatrixContainerC2C::SetToZero
void SetToZero(IndicesVectorType &data_vector)
Definition: edge_data_c2c.h:1071
Kratos::MatrixContainerC2C::GetTuplePointer
CSR_Tuple * GetTuplePointer(unsigned int NodeI, unsigned int NeighbourJ)
Definition: edge_data_c2c.h:741
Kratos::MatrixContainerC2C::FillScalarFromDatabase
void FillScalarFromDatabase(Variable< double > &rVariable, ValuesVectorType &rDestination, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:871
Kratos::MatrixContainerC2C::IndicesVectorType
vector< unsigned int > IndicesVectorType
Definition: edge_data_c2c.h:408
Kratos::MatrixContainerC2C::AllocateAndSetToZero
void AllocateAndSetToZero(ValuesVectorType &data_vector, int size)
Definition: edge_data_c2c.h:1046
Kratos::MatrixContainerC2C::GetCSRIndex
unsigned int GetCSRIndex(unsigned int NodeI, unsigned int NeighbourJ)
Definition: edge_data_c2c.h:721
Kratos::MatrixContainerC2C::MatrixContainerC2C
MatrixContainerC2C()
Definition: edge_data_c2c.h:415
Kratos::MatrixContainerC2C::AssignVectorToVector
void AssignVectorToVector(const CalcVectorType &origin, CalcVectorType &destination)
Definition: edge_data_c2c.h:1102
Kratos::MatrixContainerC2C::BuildCSRData
void BuildCSRData(ModelPart &model_part)
Definition: edge_data_c2c.h:584
Kratos::MatrixContainerC2C::SetToZero
void SetToZero(ValuesVectorType &data_vector)
Definition: edge_data_c2c.h:1091
Kratos::MatrixContainerC2C::GetInvertedMass
ValuesVectorType & GetInvertedMass()
Definition: edge_data_c2c.h:446
Kratos::MatrixContainerC2C::WriteVectorToDatabase
void WriteVectorToDatabase(Variable< array_1d< double, 3 >> &rVariable, CalcVectorType &rOrigin, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:926
Kratos::MatrixContainerC2C::AllocateAndSetToZero
void AllocateAndSetToZero(CalcVectorType &data_vector, int size)
Definition: edge_data_c2c.h:1034
Kratos::MatrixContainerC2C::Clear
void Clear()
Definition: edge_data_c2c.h:761
Kratos::MatrixContainerC2C::SetToZero
void SetToZero(CalcVectorType &data_vector)
Definition: edge_data_c2c.h:1080
Kratos::MatrixContainerC2C::WriteScalarToDatabase
void WriteScalarToDatabase(Variable< double > &rVariable, ValuesVectorType &rOrigin, ModelPart::NodesContainerType &rNodes)
Definition: edge_data_c2c.h:954
Kratos::MatrixContainerC2C::ConstructCSRVector
void ConstructCSRVector(ModelPart &model_part)
Definition: edge_data_c2c.h:464
Kratos::ModelPart
This class aims to manage meshes for multi-physics simulations.
Definition: model_part.h:77
Kratos::ModelPart::NodesContainerType
MeshType::NodesContainerType NodesContainerType
Nodes container. Which is a vector set of nodes with their Id's as key.
Definition: model_part.h:128
Kratos::OpenMPUtils::DivideInPartitions
static void DivideInPartitions(const int NumTerms, const int NumThreads, PartitionVector &Partitions)
Divide an array of length NumTerms between NumThreads threads.
Definition: openmp_utils.h:158
Kratos::OpenMPUtils::ThisThread
static int ThisThread()
Wrapper for omp_get_thread_num().
Definition: openmp_utils.h:108
Kratos::ParallelUtilities::GetNumThreads
static int GetNumThreads()
Returns the current number of threads.
Definition: parallel_utilities.cpp:34
Kratos::Variable
Variable class contains all information needed to store and retrive data from a data container.
Definition: variable.h:63
Kratos::array_1d< double, TDim >
KRATOS_CATCH
#define KRATOS_CATCH(MoreInfo)
Definition: define.h:110
KRATOS_WATCH
#define KRATOS_WATCH(variable)
Definition: define.h:806
KRATOS_TRY
#define KRATOS_TRY
Definition: define.h:109
free_surface_application.h
model_part.h
Kratos
REF: G. R. Cowper, GAUSSIAN QUADRATURE FORMULAS FOR TRIANGLES.
Definition: mesh_condition.cpp:21
Kratos::ZeroVector
KratosZeroVector< double > ZeroVector
Definition: amatrix_interface.h:561
Kratos::ZeroMatrix
KratosZeroMatrix< double > ZeroMatrix
Definition: amatrix_interface.h:559
Kratos::noalias
T & noalias(T &TheMatrix)
Definition: amatrix_interface.h:484
face_heat.model_part
model_part
Definition: face_heat.py:14
generate_convection_diffusion_explicit_element.tau
tau
Definition: generate_convection_diffusion_explicit_element.py:115
generate_total_lagrangian_mixed_volumetric_strain_element.n_nodes
int n_nodes
Definition: generate_total_lagrangian_mixed_volumetric_strain_element.py:15
ode_solve.L
int L
Definition: ode_solve.py:390
quadrature.k
int k
Definition: quadrature.py:595
rotating_cone.temp
float temp
Definition: rotating_cone.py:85
script_THERMAL_CORRECT.counter
int counter
Definition: script_THERMAL_CORRECT.py:218
sensitivityMatrix.N
N
Definition: sensitivityMatrix.py:29
sp_statistics.volume
volume
Definition: sp_statistics.py:15
tensors::i
integer i
Definition: TensorModule.f:17
openmp_utils.h
Generated on Sat Feb 3 2024 02:01:53 for KratosMultiphysics by doxygen 1.9.1