16 #ifndef __CELL_FORCES_HPP__
17 #define __CELL_FORCES_HPP__
37 boost::shared_ptr<VectorDouble>
dispX;
38 boost::shared_ptr<VectorDouble>
dispY;
59 :
MoFEM::OpCalculateScalarFieldValues(
"DISP_X", common_data.dispX) {}
67 :
MoFEM::OpCalculateScalarFieldValues(
"DISP_Y", common_data.dispY) {}
95 int nb_row = row_data.getIndices().size();
96 int nb_col = col_data.getIndices().size();
103 const VectorInt &row_indices = row_data.getIndices();
104 const VectorInt &col_indices = col_data.getIndices();
106 const int nb_gauss_pts = row_data.getN().size1();
108 C.resize(nb_row, nb_col,
false);
110 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
113 for (
int rr = 0; rr != nb_row; rr++) {
114 const double diff_nr_x = row_data.getDiffN()(gg, 2 * rr + 0);
115 const double diff_nr_y = row_data.getDiffN()(gg, 2 * rr + 1);
116 for (
int cc = 0; cc != nb_col; cc++) {
117 const double diff_nc_x = col_data.getDiffN()(gg, 2 * cc + 0);
118 const double diff_nc_y = col_data.getDiffN()(gg, 2 * cc + 1);
119 C(rr, cc) += val * (diff_nr_x * diff_nc_x + diff_nr_y * diff_nc_y);
125 &*col_indices.begin(), &*
C.data().begin(), ADD_VALUES);
161 int nb_row = row_data.getIndices().size();
162 int nb_col = col_data.getIndices().size();
169 const VectorInt &row_indices = row_data.getIndices();
170 const VectorInt &col_indices = col_data.getIndices();
172 const int nb_gauss_pts = row_data.getN().size1();
174 C.resize(nb_row, nb_col,
false);
176 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
179 for (
int rr = 0; rr != nb_row; rr++) {
180 const double diff_nr_x = row_data.getDiffN()(gg, 2 * rr + 0);
181 const double diff_nr_y = row_data.getDiffN()(gg, 2 * rr + 1);
182 for (
int cc = 0; cc != nb_col / 3; cc++) {
183 const double nc = col_data.getN(gg)[cc];
184 C(rr, 3 * cc + 0) -= val * diff_nr_x * nc;
185 C(rr, 3 * cc + 1) -= val * diff_nr_y * nc;
191 &*col_indices.begin(), &*
C.data().begin(), ADD_VALUES);
223 int nb_row = row_data.getIndices().size();
224 int nb_col = col_data.getIndices().size();
231 const VectorInt &row_indices = row_data.getIndices();
232 const VectorInt &col_indices = col_data.getIndices();
234 const int nb_gauss_pts = row_data.getN().size1();
236 C.resize(nb_row, nb_col,
false);
238 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
241 for (
int rr = 0; rr != nb_row / 3; rr++) {
242 const double nr = row_data.getN(gg)[rr];
243 for (
int cc = 0; cc != nb_col / 3; cc++) {
244 const double nc = col_data.getN(gg)[cc];
245 const double c = val * nr * nc;
246 C(3 * rr + 0, 3 * cc + 0) +=
c;
247 C(3 * rr + 1, 3 * cc + 1) +=
c;
253 &*col_indices.begin(), &*
C.data().begin(), ADD_VALUES);
286 int nb_row = data.getIndices().size();
293 const int nb_gauss_pts = data.getN().size1();
295 nF.resize(nb_row,
false);
297 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
299 for (
int rr = 0; rr != nb_row / 3; rr++) {
300 const double nr = val * data.getN(gg)[rr];
301 nF[3 * rr + 0] += nr * disp_x[gg];
302 nF[3 * rr + 1] += nr * disp_y[gg];
338 int nb_row = row_data.getIndices().size();
339 int nb_col = col_data.getIndices().size();
346 const VectorInt &row_indices = row_data.getIndices();
347 const VectorInt &col_indices = col_data.getIndices();
349 const int nb_gauss_pts = row_data.getN().size1();
351 auto t_diff_base_row = row_data.getFTensor2DiffN<3, 2>();
352 auto t_base_row = row_data.getFTensor1N<3>();
359 C.resize(nb_row, nb_col,
false);
361 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
364 for (
int rr = 0; rr != nb_row; rr++) {
365 const double curl_row = t_diff_base_row(1, 0) - t_diff_base_row(0, 1);
366 auto t_diff_base_col = col_data.getFTensor2DiffN<3, 2>(gg, 0);
367 auto t_base_col = col_data.getFTensor1N<3>(gg, 0);
368 for (
int cc = 0; cc != nb_col; cc++) {
369 const double curl_col = t_diff_base_col(1, 0) - t_diff_base_col(0, 1);
370 const double diag = val *
epsRho * t_base_row(
i) * t_base_col(
i);
372 epsL > 0 ? val * pow(
epsL, -1) * curl_row * curl_col : 0;
373 C(rr, cc) += diag + curl;
383 &*col_indices.begin(), &*
C.data().begin(), ADD_VALUES);
413 int nb_row = row_data.getIndices().size();
414 int nb_col = col_data.getIndices().size();
421 const VectorInt &row_indices = row_data.getIndices();
422 const VectorInt &col_indices = col_data.getIndices();
424 const int nb_gauss_pts = row_data.getN().size1();
426 auto t_base_row = row_data.getFTensor1N<3>();
428 C.resize(nb_row, nb_col,
false);
430 for (
int gg = 0; gg != nb_gauss_pts; gg++) {
433 for (
int rr = 0; rr != nb_row; rr++) {
435 auto t_base_col = col_data.getFTensor0N(gg, 0);
436 for (
int cc = 0; cc != nb_col / 3; cc++) {
437 C(rr, 3 * cc + 0) -= val * t_base_row(0) * t_base_col;
438 C(rr, 3 * cc + 1) -= val * t_base_row(1) * t_base_col;
446 &*col_indices.begin(), &*
C.data().begin(), ADD_VALUES);
479 if (data.getFieldData().empty())
482 const int nb_gauss_pts = data.getN().size1();
483 const int nb_dofs = data.getFieldData().size();
495 if (
type == MBVERTEX) {
496 mat->resize(nb_gauss_pts, 2,
false);
498 vec->resize(nb_gauss_pts,
false);
500 grad->resize(0, 0,
false);
504 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
505 for (
int rr = 0; rr != nb_dofs; rr++) {
506 const double diff_base_x = data.getDiffN()(gg, 2 * rr + 0);
507 const double diff_base_y = data.getDiffN()(gg, 2 * rr + 1);
508 const double val = data.getFieldData()[rr];
509 (*mat)(gg, 0) += diff_base_x * val;
510 (*mat)(gg, 1) += diff_base_y * val;
511 const double base = data.getN()(gg, rr);
512 (*vec)[gg] += base * val;
547 if (data.getFieldData().empty())
550 const int nb_gauss_pts = data.getN().size1();
551 const int nb_dofs = data.getFieldData().size();
563 if (
type == MBEDGE && side == 0) {
564 mat->resize(nb_gauss_pts, 2,
false);
566 grad->resize(nb_gauss_pts, 4,
false);
568 vec->resize(0,
false);
572 auto t_base = data.getFTensor1N<3>();
573 auto t_diff_base_fun = data.getFTensor2DiffN<3, 2>();
575 for (
unsigned int gg = 0; gg < nb_gauss_pts; gg++) {
576 auto t_data = data.getFTensor0FieldData();
577 for (
int rr = 0; rr != nb_dofs; rr++) {
578 const double base_x = t_base(0);
579 const double base_y = t_base(1);
580 (*mat)(gg, 0) += base_x * t_data;
581 (*mat)(gg, 1) += base_y * t_data;
582 for (
int i = 0;
i != 2;
i++) {
583 for (
int j = 0;
j != 2;
j++) {
584 (*grad)(gg, 2 *
i +
j) += t_diff_base_fun(
i,
j) * t_data;
612 std::vector<EntityHandle> &map_gauss_pts,
626 if (
type != MBVERTEX)
629 double def_val[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
630 Tag th_force, th_force_potential, th_grad_force;
632 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
637 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
642 MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
645 int nb_gauss_pts = data.getN().size1();
646 if (
mapGaussPts.size() != (
unsigned int)nb_gauss_pts) {
650 for (
int gg = 0; gg < nb_gauss_pts; gg++) {
681 #endif // __CELL_FORCES_HPP__