Flow123d  master-f44eb46
mortar_assembly.cc
Go to the documentation of this file.
1 /*
2  * mortar_assembly.cc
3  *
4  * Created on: Feb 22, 2017
5  * Author: jb
6  */
7 
8 #include "flow/mortar_assembly.hh"
10 #include "la/linsys.hh"
11 #include "mesh/accessors.hh"
12 #include "fem/dh_cell_accessor.hh"
14 #include <armadillo>
15 
17 : MortarAssemblyBase(eq_fields, eq_data),
18  tensor_average_(16),
19  col_average_(4),
20  quadrature_(*(eq_data->mesh))
21 {
22  isec_data_list.reserve(30);
23 
24 
25  for(uint row_dim=0; row_dim<4; row_dim ++)
26  for(uint col_dim=0; col_dim<4; col_dim ++) {
27  arma::vec row_avg = arma::vec(row_dim+1);
28  row_avg.fill(1.0 / (row_dim+1));
29  arma::vec col_avg = arma::vec(col_dim+1);
30  col_avg.fill(1.0 / (col_dim+1));
31  arma::mat avg = row_avg * col_avg.t(); // tensor product
32  tensor_average(row_dim, col_dim) = avg;
33  col_average_[row_dim] = row_avg;
34  }
35 }
36 
37 
38 void P0_CouplingAssembler::pressure_diff(const DHCellAccessor& dh_cell, double delta) {
39  isec_data_list.push_back(IsecData());
40  IsecData &i_data = isec_data_list.back();
41 
42  ElementAccessor<3> ele = dh_cell.elm();
43  const uint nsides = ele->n_sides();
44  const uint ndofs = dh_cell.n_dofs();
45 
46  i_data.dim= ele.dim();
47  i_data.delta = delta;
48  i_data.dofs.resize(nsides);
49  i_data.vel_dofs.resize(nsides);
50  i_data.dirichlet_dofs.resize(nsides);
51  i_data.dirichlet_sol.resize(nsides);
52  i_data.n_dirichlet=0;
53  i_data.ele_z_coord_=ele.centre()[2];
54 
55  for(unsigned int i_side=0; i_side < nsides; i_side++ ) {
56  // TODO: replace with DHCell getter when available for FESystem component
57  i_data.dofs[i_side] = dh_cell.get_loc_dof_indices()[(ndofs+1)/2+i_side]; //edge dof
58  i_data.vel_dofs[i_side] = dh_cell.get_loc_dof_indices()[i_side]; // side dof
59  //i_data.z_sides[i_side]=ele.side(i_side)->centre()[2];
60  Side side = *dh_cell.elm().side(i_side);
61  if (side.is_boundary()) {
62  Boundary bcd = side.cond();
64  auto type = (DarcyMH::EqFields::BC_Type)eq_fields_->bc_type.value(b_ele.centre(), b_ele);
65  //DebugOut().fmt("bcd id: {} sidx: {} type: {}\n", ele->id(), i_side, type);
66  if (type == DarcyMH::EqFields::dirichlet) {
67  //DebugOut().fmt("Dirichlet: {}\n", ele->index());
68  double bc_pressure = eq_fields_->bc_pressure.value(b_ele.centre(), b_ele);
69  i_data.dirichlet_dofs[i_data.n_dirichlet] = i_side;
70  i_data.dirichlet_sol[i_data.n_dirichlet] = bc_pressure;
71  i_data.n_dirichlet++;
72  }
73  }
74  }
75 
76 }
77 
78 
79 
80 
81 
82 /**
83  * Works well but there is large error next to the boundary.
84  */
86 {
87 
88 
89  // on the intersection element we consider
90  // * intersection dofs for master and slave
91  // those are dofs of the space into which we interpolate
92  // base functions from individual master and slave elements
93  // For the master dofs both are usualy eqivalent.
94  // * original dofs - for master same as intersection dofs, for slave
95  // all dofs of slave elements
96 
97  // form list of intersection dofs, in this case pressures in barycenters
98  // but we do not use those form MH system in order to allow second schur somplement. We rather map these
99  // dofs to pressure traces, i.e. we use averages of traces as barycentric values.
100 
101  /**
102  * TODO:
103  * - pass through the master and all slaves and collect global dofs , bcd, solution.
104  * I.e. call Nx pressure_diff not NxNx.
105  *
106  * - Is it safe to have duplicate rows in local_system?
107  * - Is it better to have more smaller local system then single big one?
108  * - use one big or more smaller local systems to set.
109  */
110 
111  ElementAccessor<3> ele = dh_cell_master.elm();
112 
113  if (ele.dim() > 2) return; // supported only for 1D and 2D master elements
114  auto &isec_list = mixed_mesh_.element_intersections_[ele.idx()];
115  if (isec_list.size() == 0) return; // skip empty masters
116 
117  //slave_ac_.setup(master_ac);
118 
119  arma::vec3 ele_centre = ele.centre();
120  double m_sigma = eq_fields_->sigma.value( ele_centre, ele);
121  double m_conductivity = eq_fields_->conductivity.value( ele_centre, ele);
122  double m_crossection = eq_fields_->cross_section.value( ele_centre, ele );
123 
124  double master_sigma = 2*m_sigma*m_conductivity *
125  2/ m_crossection;
126  /**
127  * ?? How to deal with anisotropy ??
128  * 3d-2d : compute nv of 2d triangle
129  * 2d-2d : interpret as 2d-1d-2d, should be symmetric master-slave
130  * 2d-1d : nv is tangent to 2d and normal to 1d
131  arma::dot(eq_fields_->anisotropy.value( ele_centre, ele->element_accessor())*nv, nv)
132  */
133 
134 
135 
136 
137 
138 
139 
140  isec_data_list.clear();
141  double cs_sqr_avg = 0.0;
142  double isec_sum = 0.0;
143  unsigned int slave_dim = 0;
144  uint i = 0;
145  for(; i < isec_list.size(); ++i) {
146  bool non_zero = quadrature_.reinit(isec_list[i].second);
147  DHCellAccessor dh_cell_slave(this->eq_data_->dh_.get(), quadrature_.slave_idx());
148  ElementAccessor<3> ele_slave = dh_cell_slave.elm();
149  slave_dim = ele_slave.dim();
150  if (slave_dim == ele.dim()) break;
151  if (! non_zero) continue; // skip quadratures close to zero
152 
153  double cs = eq_fields_->cross_section.value(ele_slave.centre(), ele_slave);
154  double isec_measure = quadrature_.measure();
155  //DebugOut() << print_var(cs) << print_var(isec_measure);
156  cs_sqr_avg += cs*cs*isec_measure;
157  isec_sum += isec_measure;
158  //DebugOut().fmt("Assembly23: {} {} {} ", ele.idx(), ele_slave->id(), isec_measure);
159  pressure_diff(dh_cell_slave, isec_measure);
160  }
161  if ( ! (slave_dim == 2 && ele.dim() ==2 ) ) {
162  if( fabs(isec_sum - ele.measure()) > 1E-5) {
163  string out;
164  for(auto & isec : isec_list) {
165  DHCellAccessor dh_cell_slave(this->eq_data_->dh_.get(), isec.second->bulk_ele_idx());
166  out += fmt::format(" {}", dh_cell_slave.elm().idx());
167  }
168 
169  double diff = (isec_sum - ele.measure())/ele.measure();
170  WarningOut() << print_var(diff)
171  << print_var(ele.idx())
172  << endl
173  << out;
174 
175  }
176  }
177  pressure_diff(dh_cell_master, -isec_sum);
178 
179  //DebugOut().fmt( "cs2: {} d0: {}", cs_sqr_avg, delta_0);
180  master_sigma = master_sigma * (cs_sqr_avg / isec_sum)
181  / isec_sum;
182 
183 
184  add_to_linsys(master_sigma);
185 
186 
187  // 2d-2d
188  //DebugOut() << "2d-2d";
189  if (i < isec_list.size()) {
190  isec_data_list.clear();
191  isec_sum = 0.0;
192  for(; i < isec_list.size(); ++i) {
193  quadrature_.reinit(isec_list[i].second);
194  DHCellAccessor dh_cell_slave(this->eq_data_->dh_.get(), quadrature_.slave_idx());
195  double isec_measure = quadrature_.measure();
196  isec_sum += isec_measure;
197  //DebugOut().fmt("Assembly22: {} {} {}", ele.idx(), dh_cell_slave.elm().idx(), isec_measure);
198  pressure_diff(dh_cell_slave, isec_measure);
199  }
200  pressure_diff(dh_cell_master, -isec_sum);
201 
202  master_sigma = 100 * m_conductivity/ m_crossection / isec_sum;
203 
204  add_to_linsys(master_sigma);
205  }
206 }
207 
209 {
210  // rows
211  for(IsecData &row_ele : isec_data_list) {
212  //columns
213  for(IsecData &col_ele : isec_data_list) {
214 
215 
216  double s = -scale * row_ele.delta * col_ele.delta;
217  //DebugOut().fmt("s: {} {} {} {}", s, scale, row_ele.delta, col_ele.delta);
218  product_ = s * tensor_average(row_ele.dim, col_ele.dim);
219 
220  loc_system_.reset(row_ele.dofs, col_ele.dofs);
221 
222  for(uint i=0; i< row_ele.n_dirichlet; i++)
223  loc_system_.set_solution_row(row_ele.dirichlet_dofs[i], row_ele.dirichlet_sol[i], -1.0);
224  for(uint i=0; i< col_ele.n_dirichlet; i++) loc_system_.set_solution_col(col_ele.dirichlet_dofs[i], col_ele.dirichlet_sol[i]);
225  //ASSERT_PERMANENT( arma::norm(product_,2) == 0.0 );
227  // Must have z-coords for every side, can not use averaging vector
228  loc_system_.set_rhs( -s * col_average_[row_ele.dim] * col_ele.ele_z_coord_ );
230 
231  if (fix_velocity_flag) {
232  this->fix_velocity_local(row_ele, col_ele);
233  } else {
235  eq_data_->lin_sys->set_local_system(loc_system_, eq_data_->dh_->get_local_to_global_map());
236  }
237  }
238  }
239 }
240 
241 
242 void P0_CouplingAssembler::fix_velocity_local(const IsecData &row_ele, const IsecData & col_ele)
243 {
244 
245  uint n_rows = row_ele.vel_dofs.n_rows;
246  uint n_cols = col_ele.dofs.n_rows;
247  arma::vec pressure(n_cols);
248  arma::vec add_velocity(n_rows);
249 
250  for(uint icol=0; icol < n_cols; icol++ ) pressure[icol] = eq_data_->full_solution.get(col_ele.dofs[icol]);
251  add_velocity = loc_system_.get_matrix() * pressure - loc_system_.get_rhs();
252  //DebugOut() << "fix_velocity\n" << pressure << add_velocity;
253  for(uint irow=0; irow < n_rows; irow++ ) eq_data_->full_solution.add( row_ele.vel_dofs[irow], add_velocity[irow] );
254 }
255 
256 void P1_CouplingAssembler::add_sides(const DHCellAccessor& dh_cell, unsigned int shift, vector<int> &dofs, vector<double> &dirichlet)
257 {
258  ElementAccessor<3> ele = dh_cell.elm();
259  const uint ndofs = dh_cell.n_dofs();
260 
261  for(unsigned int i_side=0; i_side < ele->n_sides(); i_side++ ) {
262  // TODO: replace with DHCell getter when available for FESystem component
263  dofs[shift+i_side] = dh_cell.get_loc_dof_indices()[(ndofs+1)/2+i_side]; //edge dof
264 
265  Side side = *dh_cell.elm().side(i_side);
266  if (side.is_boundary()) {
267  Boundary bcd = side.cond();
268  ElementAccessor<3> b_ele = bcd.element_accessor();
269  auto type = (DarcyMH::EqFields::BC_Type)eq_fields_->bc_type.value(b_ele.centre(), b_ele);
270  //DebugOut().fmt("bcd id: {} sidx: {} type: {}\n", ele->id(), i_side, type);
271  if (type == DarcyMH::EqFields::dirichlet) {
272  //DebugOut().fmt("Dirichlet: {}\n", ele->index());
273  dofs[shift + i_side] = -dofs[shift + i_side];
274  double bc_pressure = eq_fields_->bc_pressure.value(b_ele.centre(), b_ele);
275  dirichlet[shift + i_side] = bc_pressure;
276  }
277  }
278  }
279 }
280 
281 
282 /**
283  * P1 connection of different dimensions
284  *
285  * - 20.11. 2014 - very poor convergence, big error in pressure even at internal part of the fracture
286  */
287 
289 /*
290  const IsecList &master_list = master_list_[ele_ac.ele_global_idx()];
291  if (master_list.size() == 0) return; // skip empty masters
292  double master_sigma=eq_fields_->sigma.value( ele_ac.element_accessor()->centre(), ele_ac.element_accessor());
293 
294  // set mater sides
295  add_sides(ele_ac, 3, dofs, dirichlet);
296 
297  for(uint i = 0; i < master_list.size(); ++i) {
298  const IntersectionQuadrature &intersec = intersections_[ master_list[i] ];
299  const Element * slave = intersec.slave_iter();
300  ele_ac.reinit(intersec.slave_iter()->index());
301  add_sides(ele_ac, 0, dofs, dirichlet);
302 */
303 
304 /*
305  * Local coordinates on 1D
306  * t0
307  * node 0: 0.0
308  * node 1: 1.0
309  *
310  * base fce points
311  * t0 = 0.0 on side 0 node 0
312  * t0 = 1.0 on side 1 node 1
313  *
314  * Local coordinates on 2D
315  * t0 t1
316  * node 0: 0.0 0.0
317  * node 1: 1.0 0.0
318  * node 2: 0.0 1.0
319  *
320  * base fce points
321  * t0=0.5, t1=0.0 on side 0 nodes (0,1)
322  * t0=0.5, t1=0.5 on side 1 nodes (1,2)
323  * t0=0.0, t1=0.5 on side 2 nodes (2,0)
324  */
325 
326 /*
327 
328  arma::vec point_Y(1);
329  point_Y.fill(1.0);
330  arma::vec point_2D_Y(intersec.map_to_slave(point_Y)); // local coordinates of Y on slave (1, t0, t1)
331  arma::vec point_1D_Y(intersec.map_to_master(point_Y)); // local coordinates of Y on master (1, t0)
332 
333  arma::vec point_X(1);
334  point_X.fill(0.0);
335  arma::vec point_2D_X(intersec.map_to_slave(point_X)); // local coordinates of X on slave (1, t0, t1)
336  arma::vec point_1D_X(intersec.map_to_master(point_X)); // local coordinates of X on master (1, t0)
337 
338  arma::mat base_2D(3, 3);
339  // basis functions are numbered as sides
340  // TODO:
341  // Use RT finite element to evaluate following matrices.
342 
343  // Ravirat - Thomas base functions evaluated in points (0,0), (1,0), (0,1)
344  // 2D RT_i(t0, t1) = a0 + a1*t0 + a2*t1
345  // a0 a1 a2
346  base_2D << 1.0 << 0.0 << -2.0 << arma::endr // RT for side 0
347  << 1.0 << -2.0 << 0.0 << arma::endr // RT for side 1
348  << -1.0 << 2.0 << 2.0 << arma::endr;// RT for side 2
349 
350 
351  arma::mat base_1D(2, 2);
352  // Ravirat - Thomas base functions evaluated in points (0,0), (1,0), (0,1)
353  // 1D RT_i(t0) = a0 + a1 * t0
354  // a0 a1
355  base_1D << 1.0 << -1.0 << arma::endr // RT for side 0,
356  << 0.0 << 1.0 << arma::endr; // RT for side 1,
357 
358 
359 
360  // Consider both 2D and 1D value are defined for the test function
361  // related to the each of 5 DOFs involved in the intersection.
362  // One of these values is always zero.
363  // Compute difference of the 2D and 1D value for every DOF.
364  // Compute value of this difference in both endpoints X,Y of the intersection.
365 
366  arma::vec difference_in_Y(5);
367  arma::vec difference_in_X(5);
368  // slave sides 0,1,2
369  difference_in_Y.subvec(0, 2) = -base_2D * point_2D_Y;
370  difference_in_X.subvec(0, 2) = -base_2D * point_2D_X;
371  // master sides 3,4
372  difference_in_Y.subvec(3, 4) = base_1D * point_1D_Y;
373  difference_in_X.subvec(3, 4) = base_1D * point_1D_X;
374 
375  // applying the Simpson's rule
376  // to the product of two linear functions f, g we get
377  // (b-a)/6 * ( 3*f(a)*g(a) + 3*f(b)*g(b) + 2*f(a)*g(b) + 2*f(b)*g(a) )
378  arma::mat A(5, 5);
379  for (int i = 0; i < 5; ++i) {
380  for (int j = 0; j < 5; ++j) {
381  A(i, j) = -master_sigma * intersec.intersection_true_size() *
382  ( difference_in_Y[i] * difference_in_Y[j]
383  + difference_in_Y[i] * difference_in_X[j]/2
384  + difference_in_X[i] * difference_in_Y[j]/2
385  + difference_in_X[i] * difference_in_X[j]
386  ) * (1.0 / 3.0);
387 
388  }
389  }
390  auto dofs_cp=dofs;
391  eq_data_->lin_sys->set_values( dofs_cp, dofs_cp, A, rhs, dirichlet, dirichlet);
392 
393  }*/
394 }
395 
396 
397 
ElementAccessor< 3 > element_accessor()
Cell accessor allow iterate over DOF handler cells.
const ElementAccessor< 3 > elm() const
Return ElementAccessor to element of loc_ele_idx_.
unsigned int n_dofs() const
Return number of dofs on given cell.
LocDofVec get_loc_dof_indices() const
Returns the local indices of dofs associated to the cell on the local process.
double measure() const
Computes the measure of the element.
SideIter side(const unsigned int loc_index)
unsigned int idx() const
We need this method after replacing Region by RegionIdx, and movinf RegionDB instance into particular...
Definition: accessors.hh:215
unsigned int dim() const
Definition: accessors.hh:188
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
unsigned int n_sides() const
Definition: elements.h:131
bool reinit(const IntersectionLocalBase *isec)
void set_solution_col(uint loc_col, double solution)
void eliminate_solution()
const arma::mat & get_matrix()
Definition: local_system.hh:82
void set_rhs(arma::vec rhs)
const arma::vec & get_rhs()
Definition: local_system.hh:83
void set_matrix(arma::mat matrix)
void reset()
Resets the matrix, RHS, dofs to zero and clears solution settings.
Definition: local_system.cc:46
void set_solution_row(uint loc_row, double solution, double diag=0.0)
Definition: local_system.cc:98
std::vector< std::vector< ILpair > > element_intersections_
AssemblyDataPtr eq_data_
AssemblyFieldsPtr eq_fields_
LocalSystem loc_system_
MixedMeshIntersections & mixed_mesh_
void fix_velocity_local(const IsecData &row_ele, const IsecData &col_ele)
void add_to_linsys(double scale)
vector< IsecData > isec_data_list
std::vector< arma::vec > col_average_
void pressure_diff(const DHCellAccessor &dh_cell, double delta)
P0_CouplingAssembler(AssemblyFieldsPtr eq_fields, AssemblyDataPtr eq_data)
IntersectionQuadratureP0 quadrature_
arma::mat & tensor_average(unsigned int row_dim, unsigned int col_dim)
void assembly(const DHCellAccessor &dh_cell)
void assembly(const DHCellAccessor &dh_cell)
vector< double > dirichlet
void add_sides(const DHCellAccessor &dh_cell, unsigned int shift, vector< int > &dofs, vector< double > &dirichlet)
Boundary cond() const
bool is_boundary() const
Returns true for side on the boundary.
Wrappers for linear systems based on MPIAIJ and MATIS format.
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define print_var(var)
Definition: logger.hh:292
unsigned int uint
std::shared_ptr< DarcyMH::EqFields > AssemblyFieldsPtr
std::shared_ptr< DarcyMH::EqData > AssemblyDataPtr
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
std::string format(CStringRef format_str, ArgList args)
Definition: format.h:3141
#define FMT_UNUSED
Definition: posix.h:75
unsigned int dim
unsigned int n_dirichlet
arma::vec dirichlet_sol
LocDofVec dofs
double ele_z_coord_
LocDofVec vel_dofs
double delta
LocDofVec dirichlet_dofs