Skip to content

Commit

Permalink
Rebase fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller
gassmoeller committed Oct 23, 2017
1 parent 465628e commit 76a5e23
Show file tree
Hide file tree
Showing 14 changed files with 975 additions and 1,050 deletions.
48 changes: 17 additions & 31 deletions cookbooks/inner_core_convection/inner_core_assembly.cc
View file
Expand Up @@ -51,62 +51,48 @@ namespace aspect
execute (internal::Assembly::Scratch::ScratchBase<dim> &scratch_base,
internal::Assembly::CopyData::CopyDataBase<dim> &data_base) const
{
internal::Assembly::Scratch::StokesSystem<dim> *scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>* > (&scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> *data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>* > (&data_base);

Assert((scratch != NULL) && (data != NULL),
ExcInternalError());
internal::Assembly::Scratch::StokesSystem<dim> &scratch = dynamic_cast<internal::Assembly::Scratch::StokesSystem<dim>& > (scratch_base);
internal::Assembly::CopyData::StokesSystem<dim> &data = dynamic_cast<internal::Assembly::CopyData::StokesSystem<dim>& > (data_base);

const Introspection<dim> &introspection = this->introspection();
const FiniteElement<dim> &fe = this->get_fe();
const unsigned int stokes_dofs_per_cell = data->local_dof_indices.size();
const unsigned int n_q_points = scratch->face_finite_element_values.n_quadrature_points;

const typename DoFHandler<dim>::active_cell_iterator cell (&this->get_triangulation(),
scratch->face_finite_element_values.get_cell()->level(),
scratch->face_finite_element_values.get_cell()->index(),
&this->get_dof_handler());

unsigned int face_no = numbers::invalid_unsigned_int;
for (face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
if (scratch->face_finite_element_values.get_face_index() == cell->face_index(face_no))
break;
Assert(face_no != numbers::invalid_unsigned_int,ExcInternalError());
const unsigned int stokes_dofs_per_cell = data.local_dof_indices.size();
const unsigned int n_q_points = scratch.face_finite_element_values.n_quadrature_points;

//assemble force terms for the matrix for all boundary faces
if (cell->face(face_no)->at_boundary())
if (scratch.cell->face(scratch.face_number)->at_boundary())
{
scratch->face_finite_element_values.reinit (cell, face_no);
scratch.face_finite_element_values.reinit (scratch.cell, scratch.face_number);

for (unsigned int q=0; q<n_q_points; ++q)
{
const double P = dynamic_cast<const MaterialModel::InnerCore<dim>&>
(this->get_material_model()).resistance_to_phase_change
.value(scratch->material_model_inputs.position[q]);
.value(scratch.material_model_inputs.position[q]);

for (unsigned int i = 0, i_stokes = 0; i_stokes < stokes_dofs_per_cell; /*increment at end of loop*/)
{
if (introspection.is_stokes_component(fe.system_to_component_index(i).first))
{
scratch->phi_u[i_stokes] = scratch->face_finite_element_values[introspection
.extractors.velocities].value(i, q);
scratch.phi_u[i_stokes] = scratch.face_finite_element_values[introspection
.extractors.velocities].value(i, q);
++i_stokes;
}
++i;
}

const Tensor<1,dim> normal_vector = scratch->face_finite_element_values.normal_vector(q);
const double JxW = scratch->face_finite_element_values.JxW(q);
const Tensor<1,dim> normal_vector = scratch.face_finite_element_values.normal_vector(q);
const double JxW = scratch.face_finite_element_values.JxW(q);

// boundary term: P*u*n*v*n*JxW(q)
for (unsigned int i=0; i<stokes_dofs_per_cell; ++i)
for (unsigned int j=0; j<stokes_dofs_per_cell; ++j)
data->local_matrix(i,j) += P *
scratch->phi_u[i] *
normal_vector *
scratch->phi_u[j] *
normal_vector *
JxW;
data.local_matrix(i,j) += P *
scratch.phi_u[i] *
normal_vector *
scratch.phi_u[j] *
normal_vector *
JxW;
}
}
}
Expand Down
8 changes: 7 additions & 1 deletion include/aspect/free_surface.h
View file
Expand Up @@ -89,7 +89,13 @@ namespace aspect
*/
void parse_parameters (ParameterHandler &prm);

double get_theta () const;

/**
* Return the chosen stabilization term. See
* Kaus et al 2010 for details on the meaning of
* this term.
*/
double get_stabilization_term () const;

private:
/**
Expand Down
30 changes: 14 additions & 16 deletions include/aspect/melt.h
View file
Expand Up @@ -258,28 +258,26 @@ namespace aspect
public:
virtual ~MeltPressureRHSCompatibilityModification () {};

/**
* Integrate the local fluid pressure shape functions on a single cell
* for models with melt migration, so that they can later be used to do
* the pressure right-hand side compatibility modification.
*/
virtual
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
internal::Assembly::CopyData::CopyDataBase<dim> &data) const;
};

/**
* Assemble traction boundary condition terms for models with melt.
*/
template <int dim>
void
boundary_traction_melt (const SimulatorAccess<dim> &simulator_access,
const typename DoFHandler<dim>::active_cell_iterator &cell,
const unsigned int face_no,
internal::Assembly::Scratch::StokesSystem<dim> &scratch,
internal::Assembly::CopyData::StokesSystem<dim> &data);
}
/**
* Assemble traction boundary condition terms for models with melt.
*/
template <int dim>
class MeltBoundaryTraction : public MeltInterface<dim>
{
public:
virtual ~MeltBoundaryTraction () {};

virtual
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &scratch,
internal::Assembly::CopyData::CopyDataBase<dim> &data) const;
};
}


Expand Down
39 changes: 35 additions & 4 deletions include/aspect/simulator/assemblers/interface.h
View file
Expand Up @@ -55,6 +55,18 @@ namespace aspect
ScratchBase() {};

virtual ~ScratchBase () {};

/**
* Cell object on which we currently operate.
*/
typename DoFHandler<dim>::active_cell_iterator cell;

/**
* The number of the face object with respect to the current
* cell on which we operate. If we currently
* operate on a cell, this is not meaningful.
*/
unsigned face_number;
};

template <int dim>
Expand Down Expand Up @@ -90,6 +102,11 @@ namespace aspect
MaterialModel::MaterialModelInputs<dim> material_model_inputs;
MaterialModel::MaterialModelOutputs<dim> material_model_outputs;

/**
* Whether the Stokes matrix should be rebuild during this
* assembly. If the matrix does not change, assembling the right
* hand side is sufficient.
*/
const bool rebuild_stokes_matrix;
};

Expand All @@ -114,7 +131,8 @@ namespace aspect
const unsigned int stokes_dofs_per_cell,
const bool add_compaction_pressure,
const bool use_reference_profile,
const bool rebuild_stokes_matrix);
const bool rebuild_stokes_matrix,
const bool rebuild_stokes_newton_matrix);

StokesSystem (const StokesSystem<dim> &data);

Expand Down Expand Up @@ -148,6 +166,13 @@ namespace aspect
*/
std::vector<double> reference_densities;
std::vector<double> reference_densities_depth_derivative;

/**
* Whether the Newton solver Stokes matrix should be rebuild during
* this assembly. If the matrix does not change, assembling the right
* hand side is sufficient.
*/
const bool rebuild_newton_stokes_matrix;
};


Expand Down Expand Up @@ -329,10 +354,11 @@ namespace aspect
AdvectionSystem (const AdvectionSystem &data);

/**
* Local contributions to the global matrix and right hand side
* Local contributions to the global matrix
* that correspond only to the variables listed in local_dof_indices
*/
FullMatrix<double> local_matrix;

/** Local contributions to the global matrix from the face terms in the
* discontinuous Galerkin method. The vectors are of length
* GeometryInfo<dim>::max_children_per_face * GeometryInfo<dim>::faces_per_cell
Expand All @@ -346,6 +372,11 @@ namespace aspect
std::vector<FullMatrix<double> > local_matrices_int_ext;
std::vector<FullMatrix<double> > local_matrices_ext_int;
std::vector<FullMatrix<double> > local_matrices_ext_ext;

/**
* Local contributions to the right hand side
* that correspond only to the variables listed in local_dof_indices
*/
Vector<double> local_rhs;

/** Denotes which face matrices have actually been assembled in the DG field
Expand All @@ -368,7 +399,7 @@ namespace aspect
/** Indices of the degrees of freedom corresponding to the temperature
* or composition field on all possible neighboring cells. This is used
* in the discontinuous Galerkin method. The outer std::vector has
* length GeometryInfo<dim>::max_children_per_face * GeometryInfo<dim>::faces_per_cell,
* length GeometryInfo<dim>::max_children_per_face * GeometryInfo<dim>::faces_per_cell,
* and has size zero if in the continuous Galerkin case.
**/
std::vector<std::vector<types::global_dof_index> > neighbor_dof_indices;
Expand Down Expand Up @@ -402,7 +433,7 @@ namespace aspect
virtual
void
execute(internal::Assembly::Scratch::ScratchBase<dim> &,
internal::Assembly::CopyData::CopyDataBase<dim> &) const {}
internal::Assembly::CopyData::CopyDataBase<dim> &) const =0;

/**
* This function gets called if a MaterialModelOutputs is created
Expand Down

0 comments on commit 76a5e23

Please sign in to comment.