/**
* Initialize the matrices
- * #interface_up and #interface_down
+ * #interface_in and #interface_out
* used for local refinement
* with continuous
* Galerkin methods.
*/
- void initialize_interfaces(MGLevelObject<MATRIX>& interface_up,
- MGLevelObject<MATRIX>& interface_down);
+ void initialize_interfaces(MGLevelObject<MATRIX>& interface_in,
+ MGLevelObject<MATRIX>& interface_out);
/**
* Initialize the local data
* in the
* into a global matrix.
*/
void assemble(MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2);
+
+ /**
+ * Assemble a single matrix
+ * into a global matrix.
+ */
+ void assemble(MATRIX& G,
const FullMatrix<double>& M,
const std::vector<unsigned int>& i1,
const std::vector<unsigned int>& i2,
- const unsigned int level = -1);
+ const unsigned int level);
/**
* Assemble a single matrix
void assemble_transpose(MATRIX& G,
const FullMatrix<double>& M,
const std::vector<unsigned int>& i1,
- const std::vector<unsigned int>& i2,
- const unsigned int level = -1);
+ const std::vector<unsigned int>& i2);
/**
* Assemble a single matrix
* into a global matrix.
*/
- void assemble_down(MATRIX& G,
+ void assemble_in(MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2,
+ const unsigned int level = -1);
+
+ /**
+ * Assemble a single matrix
+ * into a global matrix.
+ */
+ void assemble_out(MATRIX& G,
const FullMatrix<double>& M,
const std::vector<unsigned int>& i1,
const std::vector<unsigned int>& i2,
* refinement edge, coupling
* coarse to fine.
*/
- SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_up;
+ SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_in;
/**
* The matrix used for face
* refinement edge, coupling
* fine to coarse.
*/
- SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_down;
+ SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_out;
/**
* A pointer to the object containing constraints.
*/
template <class MATRIX>
inline void
MGMatrixSimple<MATRIX>::initialize_interfaces(
- MGLevelObject<MATRIX>& up, MGLevelObject<MATRIX>& down)
+ MGLevelObject<MATRIX>& in, MGLevelObject<MATRIX>& out)
{
- interface_up = &up;
- interface_down = &down;
+ interface_in = ∈
+ interface_out = &out;
}
}
+ template <class MATRIX>
+ inline void
+ MGMatrixSimple<MATRIX>::assemble(
+ MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2)
+ {
+ AssertDimension(M.m(), i1.size());
+ AssertDimension(M.n(), i2.size());
+
+ for (unsigned int j=0; j<i1.size(); ++j)
+ for (unsigned int k=0; k<i2.size(); ++k)
+ if (std::fabs(M(j,k)) >= threshold)
+ G.add(i1[j], i2[k], M(j,k));
+ }
+
template <class MATRIX>
inline void
for (unsigned int j=0; j<i1.size(); ++j)
for (unsigned int k=0; k<i2.size(); ++k)
if (std::fabs(M(j,k)) >= threshold)
- if(mg_constrained_dofs->at_refinement_edge(level, i1[j])==false &&
- mg_constrained_dofs->at_refinement_edge(level, i2[k])==false)
- if((mg_constrained_dofs->is_boundary_index(level, i1[j])==false &&
- mg_constrained_dofs->is_boundary_index(level, i2[k])==false)
- ||
- (mg_constrained_dofs->is_boundary_index(level, i1[j])==true &&
- mg_constrained_dofs->is_boundary_index(level, i2[k])==true
- &&
- i1[j] == i2[k])
- )
- G.add(i1[j], i2[k], M(j,k));
+ if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+ {
+ if (mg_constrained_dofs->set_boundary_values())
+ {
+ if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+ !mg_constrained_dofs->is_boundary_index(level, i2[k]))
+ ||
+ (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
+ mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
+ i1[j] == i2[k]))
+ G.add(i1[j], i2[k], M(j,k));
+ }
+ else
+ G.add(i1[j], i2[k], M(j,k));
+ }
}
}
+
+ template <class MATRIX>
+ inline void
+ MGMatrixSimple<MATRIX>::assemble_transpose(
+ MATRIX& G,
+ const FullMatrix<double>& M,
+ const std::vector<unsigned int>& i1,
+ const std::vector<unsigned int>& i2)
+ {
+ AssertDimension(M.n(), i1.size());
+ AssertDimension(M.m(), i2.size());
+
+ for (unsigned int j=0; j<i1.size(); ++j)
+ for (unsigned int k=0; k<i2.size(); ++k)
+ if (std::fabs(M(k,j)) >= threshold)
+ G.add(i1[j], i2[k], M(k,j));
+ }
+
template <class MATRIX>
inline void
- MGMatrixSimple<MATRIX>::assemble_down(
+ MGMatrixSimple<MATRIX>::assemble_in(
MATRIX& G,
const FullMatrix<double>& M,
const std::vector<unsigned int>& i1,
for (unsigned int j=0; j<i1.size(); ++j)
for (unsigned int k=0; k<i2.size(); ++k)
if (std::fabs(M(j,k)) >= threshold)
- if(mg_constrained_dofs->at_refinement_edge(level, i1[j])==true &&
- mg_constrained_dofs->at_refinement_edge(level, i2[k])==false)
- if((mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==false &&
- mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==false)
- ||
- (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==true &&
- mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==true
- &&
- i1[j] == i2[k])
- )
- G.add(i1[j], i2[k], M(j,k));
+ if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+ {
+ if (mg_constrained_dofs->set_boundary_values())
+ {
+ if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+ ||
+ (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+ i1[j] == i2[k]))
+ G.add(i1[j], i2[k], M(j,k));
+ }
+ else
+ G.add(i1[j], i2[k], M(j,k));
+ }
}
}
template <class MATRIX>
inline void
- MGMatrixSimple<MATRIX>::assemble_transpose(
+ MGMatrixSimple<MATRIX>::assemble_out(
MATRIX& G,
const FullMatrix<double>& M,
const std::vector<unsigned int>& i1,
for (unsigned int j=0; j<i1.size(); ++j)
for (unsigned int k=0; k<i2.size(); ++k)
if (std::fabs(M(j,k)) >= threshold)
- if(mg_constrained_dofs->at_refinement_edge(level, i1[j])==true &&
- mg_constrained_dofs->at_refinement_edge(level, i2[k])==false)
- if((mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==false &&
- mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==false)
- ||
- (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==true &&
- mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==true
- &&
- i1[j] == i2[k])
- )
- G.add(i1[j], i2[k], M(k,j));
+ if(mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+ {
+ if (mg_constrained_dofs->set_boundary_values())
+ {
+ if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+ ||
+ (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+ i1[j] == i2[k]))
+ G.add(i1[j], i2[k], M(k,j));
+ }
+ else
+ G.add(i1[j], i2[k], M(k,j));
+ }
}
}
assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
if(mg_constrained_dofs != 0)
{
- assemble_transpose((*interface_down)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
- assemble_down((*interface_up)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
+ assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
+ assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
}
}
assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
}
+ else
+ {
+ assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
+ assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
+ assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
+ assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+ }
}
else
{