Hello!
I am trying to solve the 30-band kp model for silicon. I wrote some code
that solves this problem using Lagrangian finite elements. The picture for
the sqrt(|<psi | psi >|) wave function is in the attached file. It is clear
that there are oscillations for the ground state, which in theory should
not be there. I want to rewrite this problem for the standing wave basis,
but when I try with this basis, I do not get the correct eigenvalues and
wave functions. Tell me what to fix in my code? If this piece of code is
not enough, I can send you the full version
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/41937c29-3a3a-4424-9ae5-2c6c4d122254n%40googlegroups.com.
template <int dim>
class SineBasis : public FiniteElement<dim, dim>
{
class InternalData:FiniteElement<dim,dim>::InternalDataBase
{
public:
InternalData();
virtual ~InternalData() = default;
UpdateFlags update_each;
virtual std::size_t memory_consumption() const;
};
public:
SineBasis(const unsigned int degree)
: FiniteElement<dim, dim>(
FiniteElementData<dim>(get_dofs_per_object(degree),
1, // n_components
degree,
FiniteElementData<dim>::H1),
std::vector<bool>(1, true), // restriction_additivity
std::vector<ComponentMask>(degree, ComponentMask(std::vector<bool>(1, true)))),
degree(degree)
{
this->unit_support_points.resize(degree);
for (unsigned int i = 0; i < degree; ++i)
this->unit_support_points[i] = Point<dim>((i + 1.0) / (degree + 1.0));
// Установка флагов обновления
this->evaluation_flags = update_values | update_gradients | update_JxW_values;
// data.update_each =this->evaluation_flags;
}
virtual std::unique_ptr<FiniteElement<dim, dim>> clone() const override {
return std::make_unique<SineBasis<dim>>(*this);
}
virtual std::string get_name() const override {
return "SineBasis<" + std::to_string(dim) + ">";
}
virtual double shape_value(const unsigned int index,
const Point<dim> &point) const override
{
return basis_value(index, point);
}
virtual Tensor<1, dim> shape_grad(const unsigned int index,
const Point<dim> &point) const override
{
return basis_gradient(index, point);
}
virtual UpdateFlags requires_update_flags(const UpdateFlags flags) const override {
return flags;
}
virtual std::unique_ptr<typename FiniteElement<dim,dim>::InternalDataBase>
get_data(const UpdateFlags update_flags,
const Mapping<dim, dim> &mapping,
const Quadrature<dim> &quadrature,
internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &output_data) const override
{
auto data = std::make_unique<typename FiniteElement<dim,dim>::InternalDataBase>();
data->update_each = update_flags;
return data;
}
virtual void fill_fe_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const CellSimilarity::Similarity cell_similarity,
const Quadrature<dim> &quadrature,
const Mapping<dim, dim> &mapping,
const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
const internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
const typename FiniteElement<dim,dim>::InternalDataBase &fe_internal,
internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &output_data) const override
{
// Игнорируем ячейку и mapping, так как базис определен в эталонных координатах
(void)cell;
(void)mapping;
(void)mapping_internal;
(void)mapping_data;
for (unsigned int q = 0; q < quadrature.size(); ++q)
{
const Point<dim> &point = quadrature.point(q);
for (unsigned int i = 0; i < this->dofs_per_cell; ++i)
{
if (this->evaluation_flags & update_values)
output_data.shape_values(i, q) = basis_value(i, point);
if (this->evaluation_flags & update_gradients)
output_data.shape_gradients[i][q] = basis_gradient(i, point);
}
}
}
virtual void fill_fe_face_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const unsigned int face_no,
const Quadrature<dim-1> &quadrature,
const Mapping<dim, dim> &mapping,
const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
const internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
const typename FiniteElement<dim,dim>::InternalDataBase &fe_internal,
internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &output_data) const override
{
Assert(false, ExcNotImplemented());
}
virtual void fill_fe_subface_values(
const typename Triangulation<dim, dim>::cell_iterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const Quadrature<dim-1> &quadrature,
const Mapping<dim, dim> &mapping,
const typename Mapping<dim, dim>::InternalDataBase &mapping_internal,
const internal::FEValuesImplementation::MappingRelatedData<dim, dim> &mapping_data,
const typename FiniteElement<dim,dim>::InternalDataBase &fe_internal,
internal::FEValuesImplementation::FiniteElementRelatedData<dim, dim> &output_data) const override
{
Assert(false, ExcNotImplemented());
}
virtual bool has_support_on_face(
const unsigned int shape_index,
const unsigned int face_index) const override
{
// Все базисные функции равны нулю на границах (для условий Дирихле)
return false;
}
private:
static std::vector<unsigned int> get_dofs_per_object(const unsigned int degree) {
std::vector<unsigned int> dofs(dim+1, 0);
dofs[dim] = degree; // Все степени свободы на ячейке
return dofs;
}
double basis_value(const unsigned int index, const Point<dim> &point) const {
const double x_ref = point[0];
const double n = index + 1;
return std::sin(n * numbers::PI * x_ref);
}
Tensor<1, dim> basis_gradient(const unsigned int index, const Point<dim> &point) const {
const double x_ref = point[0];
const double n = index + 1;
Tensor<1, dim> grad;
grad[0] = n * numbers::PI * std::cos(n * numbers::PI * x_ref);
return grad;
}
const unsigned int degree;
UpdateFlags evaluation_flags;
// SineBasis<dim>::InternalData data;
};
